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A Study of a Numerical Solution to a Two- 
Dimensional Hydrodynamical Problem 


By A. Blair, N. Metropolis, J. von Neumann,* A. H. Taub & M. Tsingou 


1. Introduction. The purpose of this paper is to report the results obtained on 
Maniac I when that machine was used to solve numerically a set of difference 
equations approximating the equations of two-dimensional motion of an incom- 
pressible fluid in Eulerian coordinates. More precisely, the problem was concerned 
with the two-dimensional motion of two incompressible fluids subject only to gravi- 
tational and hydrodynamical forces which at time t = 0 were distributed as illus- 
trated in Fig. 1. 

This problem was discussed and formulated for machine computation by John 
von Neumann and others. His own original draft of a discussion of the differential 
and difference equations is given in Appendix I, and an iteration scheme for solving 
systems of linear equations is given in Appendix II. In the main body of this paper 
we shall outline the derivation of the equations employed by the computer and 
refer to these appendices for detailed discussions concerning them where necessary. 
Some of von Neumann’s difference equations were modified in the course of the 
work. The reasons for these modifications and their nature will be enlarged upon in 
the course of the discussion. 


2. The Equations of Motion and Boundary Conditions. We denote by z and y 
the Cartesian abscissa and ordinate of a point in a fixed coordinate system in a 
vertical plane oriented as in Fig. 1; that is, x and y are Eulerian coordinates. The 
velocity of the fluid at this point at time ¢ will be said to have z and y components 
u(x, y, t) and v(2, y, t), respectively. The density of the fluid will be denoted by p, 
the pressure by p, and the acceleration of gravity by g. 

The system of equations describing the motion of an incompressible fluid subject 
to the force of gravity in the verticai direction is then 


ou du 


1 dp 
(2.1) at + Uu az 


du 
Fay pas 


dv ov _ _1op 
(2.2) a + e— 


ou 
> e>- 
dy p oy 


+9 
Op 


(2.3) ae 


Op Op 
u = +1 ay 


du Ov 
A ox ehh, soem 
(24) z+ ay 0 
Equations (2.1) and (2.2) represent the conservation of momentum, equation (2.3) 
_is the incompressibility condition, and equation (2.4) states the conservation 
of mass. 
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Fic. 1.—Initial density distribution. BB, denotes the boundary separating the fluid of 
density p = 1 from that of density p = } at time / = 0. 


At any exterior boundary of the fluid the component of velocity normal to the 
boundary vanishes. That is 
(2.5) u cos (x, n) + v cos (y,n) = 0 


where cos (x, n) and cos (y, n) are the cosines of the angles between the z axis and 
the normal to the boundary and the y axis and the normal to the boundary, re- 
spectively. 

Along a curve across which there is a density discontinuity we must have the 
component of the velocity normal to the curve continuous. That is, we must have 
(2.6) [u] cos (2, n) + [v] cos (y, n) = 0 
where 

os + 4 = a 
fl=fa@',y) -f@,y) 
r + _ ~- ° ° ° ° e ° 
and x",y andz ,y represent contiguous points on opposite sides of the curve of 


discontinuities. 


3. Discontinuities. The only discontinuities present in incompressible fluid 
motion are density discontinuities and across these equation (2.6) must hold. In 
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the calculation to be described, such discontinuities are not explicitly taken into 
account. Indeed, it was one purpose of the calculation to see if this type of dis- 
continuity could be followed in time from a plot of the density contours when no 
special provision was made to provide for the discontinuities. 

Initially the density distribution assumed was that given in Fig. 1. Some pro- 
visions must be made in the difference equations to represent the spatial derivatives 
of the density on the curve BB, at time t = 0 (and at subsequent times on the 
curve into which BB, moves). We shall discuss this point subsequently. 


4. The Stream Function. It follows from equation (2.4) that there exists a 
function y called the stream function such that 


(4.1) u= —wWy, y9= 7, 

On an exterior boundary, equation (2.5) obtains and this may be written as 
(4.2) —y, cos (x, n) + ¥.cos (y,n) = 0 

If the equation of the boundary is given parametrically by 


r= 2(8), y = y(s) 


dxz\ dy a 
(z) +(@) 


then the direction cosines of the normal are 


with 


dy 


ds 


cos (2, n) 


dx 
208 n) = — 
cos (y, n) ds 
and equation (4.2) becomes 


dx dy _ 
¥e7, t Wa = 0 


That is 


y = constant 


on the exterior boundary. This constant may be chosen to be zero and the exterior 
boundary condition becomes 


(4.3) y=0 
5. Equations Determining the Stream Function and Density. Substituting from 


equation (4.1) into equations (2.1) and (2.2), we obtain 


(Vy — Wey + Vw) = Dz 
P( Wee a VWee + Vey) =<—Py + gp 
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where the subscript denotes a derivative with respect to the variable indicated. 
Differentiating the first of these with respect to y and the second with respect to z 
and adding, we have 


(5.1) (ic)e + (phw)y = ex(g — I) — py — pl = —w 
where 
(5.2) 'T = Wey — Vober 
(5.3) T= Ww — Vibe 
(5.4) "I = Wdy — Ay 
with 
(5.5) X= ber + Wy 

If we now set 
(5.6) x= — 
then 
(5.7) —(pxz)2 — (pxy)y = —o 


where w is given in terms of y and p by the right-hand side of equation (5.1), and 
on the boundary 


(5.8) x=0 


We may regard equation (5.1) and the boundary condition y = 0 as equations 
for determining y and use equation (4.1) to define u. The pressure may then be 
calculated from equations (2.1) or (2.2). 

The density may then be determined from equation (2.3), which may be written 
as 


(5.9) Pt = WyPc — Vxhy 


The system of equations (5.6) through (5.9) defines the problem to be approxi- 
mated by difference equations and to be solved numerically on the computer. 


6. The Difference Equations. For any function a(z, y, ¢) let 


(6.1) ai; = a(x,, yj, t) 
where 
x; = tAr 7=0,1,---J-1,] 
(6.2) yj = jay j=0,1,---J-1,d 
t' = hat h = 0,1, 2, --- 


and J and J are integers. Occasionally indices i + 4, 7 + 3, and h + } are used. 

The difference equations we shall consider will involve quanties y'; , xij, and 
pc; for h, i, and j given by equation (6.2). The mesh points with i = 0 or 
I(j = 0,1,---J)orj = Oor J (i = 0,1, --- J) are said to be boundary points. 
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ed. The remaining mesh points are called interior points. Equations (4.3) and (5.8) 
OZ define y,; = x2; = 0 for boundary points. Hence, we must give an algorithm for 
determining y; for interior points and p:; for interior and boundary points. This 
algorithm involves the finite difference representation of equation (5.7) for interior 
points and equation (5.9) for all points. 

Equation (5.7) is replaced by 


aay | - Pith (xi3 — xis) + pi-4.j (xij - xi-1;)] 
(6.3) 
+ Tig | Pii+4 (xi. CS xis) + pij—4 (x). _ xi ;-1)] = 4 


for interior points, that is, for 


lsisI-1 
lisjaJ-1 
where 
h h 
Qpr+4,j = pit1j + Pij 
(6.4) Qpr. +4 = pejtt + pe 
-_ and w’; is formed from y,; and p:,; as indicated in equation (A.14) of Appendix I. 
Equation (5.6) is replaced by the equation 
(6.5) vit) = wiz’ — 2atyi:; 
ons for 
| be 
h>0 
ten 
and by 
(6.6) V = Vij — Atxi.; 
‘ when 
Oxi- 
h=0 


Equation (5.9) is replaced by the equation 


h+l h *) (pi; — pias) if (y,)t5 
Pi,j eos Pij + Az 


(pes15 = p:.j) if (Wi. vs 


A A 


0) 
o| (Wy) tt? 


(6.7) : 
1,1 a (**) (pi; — pij-1) if (wei? = af (y.)M! 
1.J Ay (pij41 > pr.) if (y2) i <q 0) 

for interior points, where 

9 yes a vey h 

ed. 2(vz)i = (pz): + (Wz): 
and and a similar equation defines ( W)i rs 
) or Ona boundary such asj = J, (W,):.; = 0 from the boundary conditions, and 


ints. if initially pz = constant, it will follow from equation (6.7) that pis = Pis- 
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That is, the density discontinuity will not be able to reach the boundary j = J. 
For this reason equation (6.7) does not seem to be a suitable one for determining 
the time behavior of the density at the boundaries. It was replaced by 


ptt) ae es {pis (Wii? — — pit Wn), if (y,)r* < 0} 

(6.8) 3 Wo thy — pis odes" if (y,)i* = 0} 
at [ors e)i5! — piss (eit if (ve)%5! = 0) 

AY [hia wet Mt — ph (ye)it? if (y.)it! <0} 


for boundary points. This equation is a finite difference representation of 
pi = —(pu), — (pv)y 


just as equation (6.7) is of equation (5.9). 
Equation (6.7) differs from the finite difference form of equation (5.9) proposed 
by von Neumann, namely 


(6.9) piss = pig + At(pi ig + ~ (pu) ii 

where (pes and (pee) are evaluated from the values of p%;, vi; and Xi. j by 
substituting into the centered finite difference representation of equation (5.9) and 
the equation obtained by differentiating this equation with respect to ¢ and sub- 
stituting for p; from (5.9) and x for —y;. 

In the early calculations the finite difference form of equation (6.9) was used. 
However, it was found that near a discontinuity in the density p the values of p 
increased on the high side of the discontinuity and decreased on the low side, thus 
steadily increasing the size of the discontinuity. This unstable behavior did not 
occur when equation (6.7) was used. 


7. The solution of Equation (6.3). This equation is of the form of a set of linear 
equations which may be written as 


(7.1) Ax = 


where x is an unknown M[{= (J — 1)(J — 1)] dimensional vector, w is a known 
vector of this many dimensions, and A is a known M X M matrix. 

In Appendix II von Neumann discusses iteration schemes for solving these 
equations. He concludes that if A is a positive (or negative) definite matrix [the 
matrix A of equation (6.3) is negative definite] with a largest proper value less 
than or equal to 6 and a smallest one greater than or equal to a, then the “best” 
(in the sense defined in Appendix II) iteration scheme is given by the equation 


(7.2) n't) = Wess |, - n* ee — (w — An | + i. 
where 
by = | 
(7.3) et ; 
68-7 St 
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2a 
74) = ——— 
4 “ a+b 
In order to apply this scheme, bounds for the lowest and highest proper values, 
the numbers a and b, must be determined. 
Using the equations given by von Neumann in Section 15 of Appendix II, we 
may set 


= 6(-— ont © 4 | tet F 
a= 44 (ate sin Ti + (ay? sin i) 


=P Re sg Sa ee 
b = 4b (nts cos 51 of (ay)? °° 3) 


where @ and b are lower and upper bounds, respectively, of the density. That is 


(7.6) 0<dés5p:;56 
for all 7 and j. 


8. The Flow Diagram. A condensed flow diagram is reproduced as Fig. 2. Be- 
fore operation, initial values of y;,; and p;,; at time h = 0 are stored. The values 
of p;,; were prepared as follows: If 7, 7 labels a mesh point which does not lie on 
the density discontinuity, the value p;,; = (xi, y;). If z;, y; lies on the density 
discontinuity, we define 


pig = Zlo(2i, Yisr) + (2, Ysr)) 
y;,; was taken to be zero initially. When the routine is started, part A is traversed, 
which sets h = 0 and optionally prints out the initial values of y and p in a format 
uniform with subsequent results. 

Part B computes the values of as. the right side of equation (5.1), for all 
values of 7, 7 corresponding to interior points, as needed in equation (7.1). The 
box with the sole notation 7, 7 indicates an induction loop repeatedly using the 
program of the box to right of it for all appropriate values of 7, 7. All derivatives in 
the formula for w;,; are computed by taking the difference of the values of the 
function at lattice points on each side, for example 


1 
(vz)is = DAz (Wis.5 — Vi-r,s) 


except in certain cases near the boundary where one of these quantities does not 
exist, and then a one-sided derivative, e.g., 


1 , _ 
Az (Wigs. — Vii) for + =0,j7 = 0,1,2, --- ,J 


is used. 

Part C computes n, the first term in the sequence of vectors n’ to be constructed 
converging to x. If h = 0, then 7’ is made zero, which is a reasonable estimate in 
the case where the liquid starts moving from rest. After the motion has proceeded 
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Fic. 2.—Condensed flow diagram. 


one or more time intervals, and so h > 0, then 7’ is given the value 
(1/At)(y* — y*™) 


which is approximately x" ', a good start on a sequence to approach x’. 

Part D solves equation (6.3) by means of the iteration and mean procedure 
characterized by equations (7.2), (7.3), and (7.4). There is an 2, j induction loop 
inside a larger k loop. For each value of k the vector n'*' with components nit is 
computed by equation (7.2), and max (nit? — ni.;) is computed. If this maximum 
is above a predetermined constant 7), then k is increased by 1, b+: is computed 
from equation (7.3), and the 7, 7 induction loop is entered again to compute the 
next term in the sequence of 7° values. When a vector n’ is obtained which is 
sufficiently close to n°’, then it is considered to be x", and the control passes to 
part E. 

Part E computes the new values y'*" by equations (6.5) and (6.6), and Part F 
computes the new values p'** by equations (6.7) and (6.8). The time index h is 
then increased by 1, the results are optionally printed, and the control passes 
again to Part B. 
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The finite difference equations used in the numerical computation according to 
= the flow diagram are: 
Part A. Part A has no formulae. 
Part B. Part B is an induction loop for computing all interior points of w;,; as 
° h A 
a function of y;,; and p;,; . The formulae used are 


(B.1) 2( Ax) (We)i,j = Via. - Vis 
(B.2) 2( Ay) (Wis = Wiss — Wis 
(B.3) 4( Ax) (Ay) (Wzy)*.3 = ¥o41.541 — Vi-1541 - $e42.5-2 + ¥e-1.5-4 





The next three formulae are computed as a five-cycle induction loop for (7’, 7’) = 


(i,j), (@+ 1,9), @— 1,9), 49 + 1), (79 — 1) 





(B.4) (Ar)*(Wee)i.9 = Wiss — Wig + Vay 

(B.5) (Ay) (Wow)inae = Wirsrga — Wire + Vir je 

(B.6) (Ax) Nir, 57 = (Ar)*(ex) irs + a( Ay)? (Yay) e57 

where a = (Ax/ Ay)’. 

(B.7a) 2(Ax)*(Ne)i,j = (Ar) Migs — (Ar) Mas if i 1,1 -1 

(B.7b) (Ar)*(Az)?,; = (Ar) *Naa,5 — (Ax); iff i=1 

| (B.7c) (Ar)*(Az)?,; = (Ar)*rb,; — (Ar) A; if i=J—1 
| (B.8a) 2(Ax)*(Ay)(Ay)t.5 = (Ar) Mijn — (Ata if 7 ¥1,J-1 

(B.8b) (Ax)*(Ay)(dy)i,5 = (Ar) Ai,j41 — (Ax) "Ai, ; if j=1 
o (B.8c) (Ar)?(Ay)(Ay)?,; = (Ar)*r 5 — (Ar)Ab ja if j=J-1 

(B.9) 2(Ax)(pe)i,s = Pisa. — Pia.s 

(B.10) 2( Ay) (py)i.i = Pritt — Pisa 

8(Ax)*(Ay)'T?,; = [2( Ax) (Ye)? sJ4(Ax) (Ay) (Yen)? 5] 

(B.11) — 4[2( Ay) (Wy) i,sJ[( Ar)? (Per)? 5] 
dure 8( Ax) (Ay)? 7; = 4[2( Ar) (ve)? A(Ay)* (Yow) 2.3] 
oh (B.12) — (2(Ay) (oy)} sI4(Ax) (Ay) (Yer)? 
num 4(Ax)*( Ay) T?,; = [2(Ax) (We)? sJ2(Ax)*(Ay) dy)? i) 
= (B.13) — [2(Ay) (vy)? M2(Ax)*(Ae) 4s] 
ch is 16(Ax)*(Ay)wt,; = [2( Ax) (pe)%,s8(Ax)*(Ay)'Ti,; — 8(Ax)*(Ay)g] 
st | (Bis) + af2( Ay) (oy)¢.s18(Ax)(Ay)* *Fi.s] + (oe, sIl4(A2)*(Ay)*T,] 
rt F where 
h is 


Ax 2 
sses wh Een 
; q (*) 
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Part C. Part C is an induction loop, with respect to 7, 7 of all lattice points to 
compute the first trial value, 7;,;, at cycle h for use in the iteration process. The 
formula is 


(C.1) Ardy 9°; = | ae (viz —Vis) if h> | 
9 if h=0) 


Part D. Part D is an induction loop with respect to k, with a smaller 7, j in- 
duction loop for each k. This is the solution of the difference equation by the iterative 
and mean method. The actual formulae are 


(D.1) Qpisss = pis + pt4.5 
(D.2) Qpi-4.5 = pis + pir. 
(D.3) Qapt i414 wt a(pi,; + pr.j41) 
As before 
wi (#) 
Ay 

(D.4) 2ap:,j-+ = a(pi,; + pi.) 
(D.5) 04.5 = Qwiss.s + Wiss + Lapijey + Lapis 

—2(Ax)*(Ay)(An')i,5 = [2et4y.s]1(Ax) (Ay) ni 41,3) 
(D.6) + [2pis, J[( Ax) (Ay) nia.s) + [2api, jss][( Ax) (Ay) nia] 


+ [2api ;-s][(Ax) (Ay) n¢.ja) — o:,,[(Ax) (Ay) nf.,5] 


(az) (Ay) (Fn*)s,j = (Ax) (Ay)af 4] + 2 RB | ({-2(Ar)*(ay) An") 
2 | (Az) 
(D.7) 
- ; [16(Ar)*(Ay)o4 J) 
where d = -, ef. Eq. (7.4). 


(Qbesil[( Aa) (Ay) (Fn*) <4] — (Ax) (Ay) nia} ) 


(D.8) (Ax)(Ay)ni3 = + [(Ar)(Ay)nij] if k= 0 
(Ax)(Ay)(Fn’):,; if k= 0 ) 
(D.9) bh = 1 for k =@ 
1 - 
(D.10) Dest = 2—(i—o for k > 0 


Part E. Part E is an induction loop with respect to 7,7 and is used to compute 
vi} for all points as a function of yi; and xi; 


2At 


‘ae ae 
eee a = WSs — Taya) 


[(Ax) (Ay) x.) if h>0 
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(E.1b) 


Vii = Vig - 


At 
(Ax) (Ay) 


[((Ax) (Ay) x’.,] 


h=0 


Part F. Part F is an induction loop with respect to 7, 7 and is used to compute 
hl , ° h h+1 h 
pi; asa function of ¥i,;, Wis , and p;;. 





A 


h at { (pis — pi-as) if (Wy) < 0 | b+ 
= 5 + a h+4 t ¢ vii 
U \ (pizs, j — pi. ji) if (Wig 20 ) 
h h : A+} 
Bt At (pi,5 — pis) if (We) 2 b (y yhet 
Ay (a! if ( Mh f z)i.j 
pi. i+ — Pi i) I V2) i.3 0 | 
for «+=0,1,--- ,J 
j3=1,2,---,J—1 
wan, At | Pt WITH - pass if (WH) <0 | 
pis = Pj +} 
| 


piss, ivy) th, 4 £e: * wy) if (y, x = > 0 | 


(F.1b) : . 
_ Al pi. (ve)? ~gr pij- (We) i i,j- . if (vz ris & 0 


~ Ay | oe zes (y,) a, a pi (we) if (Wz) ye < 0 } 
for 1=1,2,---,J7-—1 
j=0,J 


9. Stability and the Choice of At. The behavior of the solutions of equation 
(6.7) with regard to stability is similar to that of the corresponding equation in one 
dimension with a constant velocity of propagation which will be taken to be positive. 
Then the equation of conservation of mass becomes 


Pc = —Upz 


which has the finite difference representation 


b+ h At h h 
i = ee Ag hes = pj-r) 
2 
or 
(9.1) py = (1 — ap; + apis 
where 
udt 
(9.2) a= 
O Be 
Equation (9.1) has solutions of the form 
(9.3) p;' = exp E 7 (jAx — aha) 
where 


exp | —'" pat | =1+ a exp| —"7" az ~ | 
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and hence 
(9.4) lex (7 a ra 1 — 4a(1 — a) cos” es 
, p 7 | = a a T oF 
Thus 8 will be real, and the equation will be stable if and only if 
asl 
That is, if 
Ats - 
u 


The value of At in the two-dimensional calculations was chosen so that 


Ax 
< 
and 
Ay 
=< = 
Ive] Sy 


The units used were such that a change in At could be accomplished by changing 
the value of the constant representing g the acceleration of gravity, and scaling y. 
It follows from equation (9.4) that up to second-order terms in Ax 


(9.5) p= ul1—é(1— ut) aol 


Hence the elementary solutions of equation (9.1) of *the form of (9.3) may be 
written as 


(9.6) pj = exp ke (jar — hatu) | exp | -i (=) Az*ha(1l — «) | 


The first factor of this expression may be written as 


exp E (x — wt) | 


with « = jAz and ¢t = hAt. This is a solution of the differential equation. Thus 
the second factor on the right-hand side of equation (9.6) shows how each of these 
elementary solutions of the differential equation is distorted when that equation is 
replaced by the finite difference equation (9.1). 

Von Neumann in a personal communication to S. Ulam (ef. Appendix III of 
reference 1) has used the term “‘pseudo-diffusion” for the distortion associated with 
an initial distribution of density 


ie 1 for x20 
1. ~ \0 for «<0 


He has shown that if this function is taken as an initial condition for the difference 
equation (9.1), then the 0 values of p advance and the 1 values of p recede to the 
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i> 


Fic. 3.—Gravity flow of two incompressible inviscid fluids at various cycles h. 


right with a velocity 


and at the same time a pseudo-diffusion-mixing region forms around the interface 
of advance whose width is essentially measured by 


ba ~ V/(1 — a)x-Ax 
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Since a was made close to one by the choice of At, the region of pseudo-diffusion 
was small for the calculations reported here. 


10. Computation Time. The results reported in Section 11 of this paper were 
run on Maniac I with J = 15 and J = 38, that is, with 624 lattice points (518 
interior points). The program required 300 words (600 orders of code exclusive of 
print routines and exclusive of orders necessary for moving information to and 
from the magnetic drum because of the limited electrostatic storage capacity of 
Maniac I). About 3750 words of dynamic storage were required. 

The time required for running one time cycle of the program on Maniac I was 
18 seconds for each iteration cycle plus 100 seconds for all the rest of the program. 
The iteration process converges so as to give accuracy in an additional decimal 
place every 10 minutes, so that one time cycle requires about an hour for six-place 
accuracy (about 200 iterations) or a half hour for three-place accuracy (about 100 
iterations). About 40 per cent of this time, however, is used in transfers to and from 
the magnetic drum, so that this much time is to be charged to the fact that a 
4000-word problem was being run on a machine with 1000-word random-access 
memory capacity. 


11. The Results. The results of the computations done on Maniac I are sum- 
marized in Figs. 3, 4, and 5, where the lattice points occur at integer values of 


5 






— 





11> 


Fic. 4.—Velocity field at h = 10. 
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— eS 





i> 


Fic. 5.—Velocity field at h = 60. 


the abscissae (7) and ordinates (j). In the first of these the locus of 
p = 4(pupper + Plower) = 0.5625 is shown for various values of h, where hAt is the 
time elapsed since the start of the calculations. Curves are shown for h = 0, 10, 20, 
30, 40, 50, 55, and 60. 

Figure 4 shows the velocity field at h = 10, as well as the loci p = 0.5625, and 
p = 0.3875, and p = 0.7375. The latter two curves are the loci at which 30 and 
70 per cent, respectively, of the initial difference in density are achieved. The band 
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covered by these two curves gives a measure of the pseudo-diffusion phenomenon. 
It is evident from Fig. 4 that the zone of pseudo-diffusion is less than one mesh 
length in thickness. 

Figure 5 is similar to Fig. 4 but in this case h = 60. Now the pseudo-diffusion 
phenomenon has increased but on the whole is still confined to a region of the order 
of two mesh lengths. 

If Ax = Ay is taken to be 1 centimeter, then the total time covered by the 
calculations (60 At) would be 0.339 second. The maximum speed attained by the 
fiuid is 0.0184 centimeter/second. 

In Fig. 5 the density distribution in the lower right-hand corner seems to have 
a somewhat anomalous behavior. This may be due to the fact that equation (6.8) 
was only applied to the boundaries 7 = 0 and j = J and not to the boundaries 
t= Oandi = 7. 


12. Concluding Remarks. The most time-consuming part of the calculations 
performed was that devoted to the computation of x; . The subsequent computa- 
tion of the velocities u,; = —(W,)?,; and vi,; = (Wz)t.;, and the density p/,; took 
relatively little time. However, the behavior of the density was most sensitive to 
the type of integration formula used. It is not yet clear that the formulas actually 
used were the best ones from the point of view of minimizing the zone of pseudo- 
diffusion. 

It is expected that the use of equation (6.8) for interior points as well as boundary 
points or other devices will keep the region of pseudo-diffusion small enough so 
that calculations on moving incompressible fluids with moving interfaces can be 
made in Eulerian coordinates. If this conjecture would prove to be correct it would 
be possible to use Eulerian coordinates and avoid the main difficulty of working in 
Lagrangian ones; namely the necessity of introducing a new Lagrangian mesh 
periodically because neighboring particles do not remain neighboring particles. 


APPENDIX I* 


The differential equations are: 


Interior: 
(1) Ue + wu, + vu, = “7 
p 
1 
(2) ve + uv, + vv, = —5 Py py, + 9, 
(3) pi + Upz + vp, = 0, 
(4) uz+v, =0 
Boundary: 
(5) cos (x, n) u + cos (y, n) v = 0. 


* Although the material in Appendix I and Appendix II was left by von Neumann in the 
form of handwritten notes and not in a form intended for wide distribution, the value of its 
content is thought to justify its inclusion in this paper. 
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From (4): 
(6) u= —yy, v= ¥:. 


(6) replaces (4). 
( Now (5) becomes: 


—cos (x, n) fy, + cos (y, n) ¥, = O, 


i.€., 
¥.:~, = cos (x, n):cos (y, n), 
i.e., 
¥.,¥, is purely normal, 
Le., 
Vien = 0 (the tangential derivative of y vanishes). 
This means that 
¥y = C (= constant) 
along the boundary. Now replacing y by ¥y — C\does not interfere with (6) (y’s 
defining relation). Hence C = 0 may be assumed, i.e.: 
p (7) ¥=0 
; (on the boundary). (7) replaces (5). 
Now only (1), (2), (3) are left. These become: 
1 
Vy + Wy Vay — ve Vy = - Pz ; 


1 
Wet = Vy Vez + vz Wey -” ~~? Py + g; 


_. aoe Wy Pr + Wz Py _ 0, 


(8) ley + I(y, Wy)] = Bee 
(9) plvic + I(¥, v2) — g] = —Py, 
(10) pr = —I(y, p). 


(8), (9), (10) replace (1), (2), (3). p can be eliminated between (8), (9), by 
forming (8), + (9).. This gives 


Lele + I(y, Why + lelez + I(¥, Wz) _ g\}z = 0, 


1.€., 
(11) (piz)s + (Wry )y _ —{—gpz + [pl (¥, Vz). +> [el (y, Vy) ly} - 


(11) replaces (8), (9) (with p eliminated). 

Thus the entire system now consists of (10), (11) (interior) and (7) (bound- 
ary), while (6) is merely a definition. 

Rewriting (10), (11) and (6): 





Interior: 





(12) Ty, = —{[e(—g + Il, vel)]e + lol (¥, Wy) I, 
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where 
Lx = (pxz)2 + (pxy)y » 
(13) 
pr = —I(y, p). 
Boundary: 
(14) y= 0. 
This suggests the following procedure: Put 
(15) x= -%. 
Then: 
Let ¥, p be known. Then y, , p; are obtained by this procedure: Put 
(16) wo = {ol—g + I(y, ve) Je + [ol (¥, Wy), - 
Determine x from this elliptic boundary value problem: 
(17) Lx = ®, 
where 


Lx = (pxz)2 + (pXxy)y 5 
and on the boundary 


(18) x = 0. 
Then: 
(19) He = —X, 
(20) pe = —I(¥,p). ~° 


The expression (16) for w may be rewritten: 


oe = pxl—g + I(¥, ¥2)] + pll(y¥, Vz) |. + pyl (y, vy) + pll(y, Wy) ly - 


Now 


[T(¥, Vz) |. - I (vz, Wz) 3 I(y, Wrz) = I(y, Wrz), 
[T(y, Vv) ly si I(Wy, vy) + I(y, Vw) = I(y, Vy). 
Hence 
(21) ao = pi—g + I(y, ¥z)] + pyl (¥y, Wy) + pl(y¥, Wer + Vwy)- 


Equation (21) replaces (16); it is more convenient for numerical calculation. 


This is a more detailed expression for w: 


(22.1) = ver + vw; 
(22.2) 'T = Vey — Ver 5 
(22.3) "Tl = vw — Vey ; 
(22.4) "I = Way — Wars, 


II 


(22.5) w = pl(g +I) +p/1 4+ pl. 











be 


fc 
fc 
fc 





NUMERICAL’ SOLUTION TO TWO-DIMENSIONAL HYDRODYNAMIC PROBLEM 163 


In addition to this the right hand side of (20) is the negative of 























(23) ‘I = I(¥, p), 
where in detail 
(24) 'T = Yeby — Vopr - 
Thus the relevant equations are these: (17) with (22.1)—(22.5) and (on the 
boundary) (18), and then (19) and (20) with (24). 
Now introduce a finite lattice for x, y, t: 
(25.1) z= 27; = 1 Az, 1=0,1,---,J-1/, 
(25.2) y = y; =f Ay, j = 0,1,---,J —1,Jd, 
(25.3) t=f =h At, h = 0,1, 2,--- 
Occasionally indices 7 + 43, 7 + 3, h + 3 are also used. For any quantity 
(26.1) a = a(z, y, t). 
The following convention is used: 
(26.2) aij = a(z;, ys, 0’). 
(17) and (22.1)-(22.5) are needed for i = 1,---, 7 — 1;j7 = 1,---, J — 1 
only. (18) is needed for the other 7, 7 only: i = 0, J;7 = 0, 1,---, J —1,J or 
i= 0,1,---,7 —1,1;7 = 0, J. (19), (20) are needed for all 7, 7:7 = 0, 1,---, 
I—1,1;j5 = 0,1,-->,J — 1, d. 
The entire system of equations can now be rewritten as follows: 
vy’; is defined for i = 1,---, J — 1;j = 1,---, J— 1. 
pu; is defined for i = 0, 1,---, 7 — 1,1;7 = 0,1,---, J — 1, J. 
: In (A.1)-(A.14) * 
i=1,---,J-—1; j=i1,---,J—1. 
(A.1) (Ve)is = Views — Vig 
‘ fort ¥ 1,J —1; 
fori = 1 omit term : 
j fori = J — 1 omit term-__. 
f (A.2) (Vy); _ Vin ow vii 
| forj = 1,J — 1; 
' for 7 = 1 omit term___; 
for j = J — 1 omit term _. 
(A.3) (Vey) ti _ Visas = Vins Ve41j—1 + ¥i-1j-1 





* The bar notation includes the required constant times Az or Ay combinations. 
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fort #1,] —1;j7 ¥1,J —1; 
for 7 = 1 omit terms___; 
forz = J — 1 omit terms 
for 7 = 1 omit terms_ __; 
for] = J — 1 omit terms.... 


In (A.4)-(A.6): 





’ 





(7, 7’) _ (2,7), (¢ + 1,9), (¢ oy l,j), (4,3 +1), 


for? = 1,] —1;j #1,J —1; 

for 2 = 1 omit term___; 

for z = I — 1 omit term___; 

for 7 = 1 omit term_ __; 

for 7 = J — 1 omit term.... 

(Henee?’ = 1,---, 7 —1;7 = 1,---,J — 1.) 














(A.A) (Yer)i's? = Virsa — Wi + Vray 
for?’ ¥ 1,7] —1; 

for 7’ = 1 omit term___; 

for 7’ = J — 1 omit term 

(A.5) (Day) i" -_ Wh ita aa yp’) 5" + Wr jr 
forj #1,J —1; 

for 7’ = 1 omit term___; 

for 7’ = J — 1 omit term __.. 

2 yA > 5 > h 
(A.6) Nir" = (Vor) i" 5" + a(Yyy) i's’ 
where 

iin Ax\’ 
Ay] 
(A.7) (Xz) 45 = Neways - Nia 


fort # 1,2 —1; 


for 7 = 1 replace the index 7 — 1 by 1 and double the entire expression; 
fort = J — 1 replace the index 7 + 1 by J — 1 and double the entire expression. 


(A.8) (Ky) 33 = Nesaa — Kijaa 
forj ¥ 1, J — 1; 


for j = 1 replace the index 7 — 1 by 1 and double the entire expression; 
for 7 = J — 1 replace the index 7 + 1 by J — 1 and double the entire expression. 


(A.9) (Be)ti = pisij — pias. 
(A.10) (By)is = Pisnt — pis. 
(A.11) Ti; _ (Ye) tj(Vey) ts —_ (Vy) ts(Wer) ti « 


(A.12) Tis = (Ve)is(Vw)is — (Wv)ts (Gav) is « 
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for 
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(A.13) Ti; = (We) ti(Xy)ts — (Vy) ti(Xe) 4; - 
(A.14) ai; = (B)is (—9 + *Th;) + a(p,)3, 7; + wh; °T;, 


where 9 = (Ar)’ Ayg [for a, ef. (A.6)]. 
In (B.1)-(B.6): 


This definition of x); [i.e., (B.6)] is, however, implicit. 


(B.1) ‘ots = Bes + Bess. 

(B.2) ts = Bay + Bi-ws. 

(B.3) “ots = a(Bts + Besa). 

(B.4) ‘si; = a( pi; + Bes-1). 

(B.5) "ais = (Cats + *ats) + aah; + ‘a4;) (for a, ef. (A.6)]. 
(B.6) Gis Xisas + Gi; Xia “Gt; Khia + ‘ot; Ka — Pst; xb; = at 





forix~1,27-—1;j3#1,J—-—1; 
for i = 1 omit the term___-; 
fori = J — 1 omit the term___; 


for 7 = 1 omit the term..__; 


for] = J — 1 omit the term.... i 
The term with the double underscore should be C x; , the corrected z°; . . 

In (C): ’ 

t=1,---,2-1; jgH=l,---,J—1. 
(C) Vis) = Viz — 4b xi; 
where 
as At _ 
Ax Ay” 


™ In (D.1-(D.3): 
é=0,1,---,7-1,7; j=0,1,---,JI—1,J. 
As in (C): 
m. hee eae 
(D.1) Az Ay 
a = $b Yijun + Vin — Vii — WIAD), 








fori ¥ 0, J andj + 0, J; 





Ses St Ball 


for i = 0, J omit the entire expression; 
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for 7 ~ 0,7 and J = 0 omit the terms and factor - 





for 7 ~ 0, J andj = J omit the terms and factor—_.. 


(D.2) 6 = $b( Wins; — With; + vias + vith) 





fori ~ 0, J andj ~ 0, J; 

for 7 = 0, J omit the entire expression; 

for 1 = 0 omit the terms and factor__; 

for 7 ~ 0, J andi = J omit the terms and factor___. 


pis = pis(1 — @ — B) + dpinis(a + 0°) + $ptaj(—a + a’) 





(D.3) + 4piier(B + B) + 4pti-a(—B6 + 8’) 








h h h h 
+ F(pisaszs — pitas — pi-rjai + pPi-1j-1) — a8 








for: ~ 0, J andj ¥ 0, J; 


for i = 0, J omit the terms : 





for 7 = 0, J omit the terms—_W. 


[Note: The points with 7 = 0, J and j = 0, J (together!) may be bypassed.] 
Alternatively, in place of (D.3): 


e=Sgna,n=Sgnée . 


h+1 h ph h | h h | 
Pig, = Pig + (Pires — pis) |a| + (pisen pis) | B | 





D3 
( ) + (Pisei+n — Pites — Prien + pis) |a||6| 








fori ~ 0, J and J ~ 0, J; 

for 7 = 0, J omit the terms___; 

for 7 = 0, J omit the terms___. 

[Note: The points with 7 = 0, J and j = 0, J (together!) may be bypassed.] 


APPENDIX II 


1. The purpose of this paper is to find a rapidly converging iterative method for 
the solution of linear equation systems, and quite particularly of those which arise 
from the difference equation treatment of partial differential equations of the el- 
liptic type [2nd order, s (= 2, 3,---) variables]. Sections 2-6 are introductory. 
The method will be described and discussed in Sections 7-13. The application to 
the (elliptic) differential equation case will be made in Sections 14-15. Some com- 
parisons will be made. The results are summarized in Sections 11, 13, 15. 

2. Consider a system of n linear equations in n variables, written vectorially: 


(1) AE =a. 





ch 
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Here @ is a known nth order vector, A a known nth order matrix, — the unknown 
nth order vector. In order that the problem be meaningful, A must be non-singu- 
lar. This will be assumed. 

An iterative method is based on a correction step, which replaces a &, that 
may not solve (1), by a &', that, in some suitable sense, should more nearly solve 
(1). This correction step should be a linear operation F applied to the two nth 
order vectors &, a, i.e., to the 2nth order vector {&, e}. It then produces the nth 
order vector ': 


(2) ef = FIE, al. 
It is convenient, to put in place of the nth order vector &' again a 2nth order vec- 
tor, namely {&', a}. In this case let us write EZ in place of F: 
(3) (', «| = ELE, al. 
Thus £ is a 2nth order matrix. 
Note that the linearity of F means that it can be written as follows: 


(4) F{E, a} = GE + Ha, 


where G, H are nth order matrices. 
(2), (3), (4) mean that the 2nth order matrix / can be written as a 2nd order 
hypermatrix of nth order matrices, as follows: 


, __ (G\H 
(5) E = (Gla). 


Here, O, J are the (nth order) zero and unit matrix, as usual. 

3. A minimum requirement to be imposed on a correction step in the sense of 2 
is this: If &* is the solution of (1), then the correction should leave — = & un- 
changed, i.e., produce &' = (= &*). This is the “weak” condition. A reasonable 
further requirement is that if — is not a solution of (1) (€' ¥ &*), then the correc- 
tion should change £, i.e., produce a &' # &. This is the “strong” condition. That 
is, the weak (strong) condition requires that & = &* be sufficient (necessary and 
sufficient) for & = &. 

By (2), (4) & = & means 


(6) (I — G)t = Ha. 
The weak condition requires, that (1) imply (6), i.e., that always 


I — G = HA, 
G = 1 — HA. 


The strong condition requires, in addition to this, that (6) imply (1), i-e., in 
view of (7), that 


HAE = Ha 


imply 


A&E= a. 
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This means obviously that 


(8a) H is non-singular. 
Equivalently: 
(8b) 0 is not a characteristic root of H. 


Since A is non-singular, non-singularity of H is equivalent to that of HA, ice., 
iby (7)] of I — G; i.e., equivalent to this: 0 is not a characteristic root of J — G, 
or equivalently : 


(8c) 1 is not a characteristic root of G. 


To begin with, we will only stipulate the weak condition, i.e., (7). 

4. The ordinary iterative procedure consists of repeating the basic step (2) 
successively, and to expect that the sequence so generated will converge to the 
solution &* of (1), irrespective of the starting point &. 

Disregarding for the moment the question of convergence, the sequence in 
question, &°, &', &,---, is defined by 


ge =, 
et — Pie* a} (k = 0, 1, 2,---).] 


j 


(9) 


In view of (2), (3), the second equation of (9) can be written 
fe" a} = Efe’, a} (k = 0, 1, 2,---) 
and hence (9) is equivalent to 
(10) fe ag = Be} - (k = 0, 1, 2,---). 
We know that for ¢ = & [* the solution of (1), ef. 3] all & = &, hence (10) 
gives 
{&*, @} = E'{e*, a}. 
Hence (10) is equivalent to what is obtained by subtracting this equation from 
it, i.e., to 
(e — &*, 0} = EE — &, 0}, 
i.e., in view of (5) to 
(11) B— & = GE — &) (k = 0, 1, 2,---). 
Note that (10) is an effective calculational procedure, while (11) is not, since 
it contains the unknown &*; however, some proofs and evaluations can be more 
advantageously based on (11). 
It is well known that frequently the convergence properties of a sequence can 
be significantly improved by replacing each element of the sequence by a suitable 


mean of itself and the preceding elements of the sequence. In this sense, one might 
replace the sequence &’, &', &’,--- by a sequence n’, n’, n’,- --, where 


k 
(12) n’ ” > anit’ (k = 0, l, 2,°°*), 
l=0 
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with a suitable set of coefficients a,,;. The characterization of the n‘* as means (of 
the ¢) makes it natural to require 


k 
(13) > aw = 1 (k = 0, 1, 2,---). 
l=0 


This condition can also be obtained from the natural requirement that for — = &*, 
when all ¢* = &, there shall also be all n* = &*. At any rate, we stipulate (13). 
The characterization of the n‘ as means might also suggest the requirement that 
all ax; 2 0, but we will not impose it; indeed, the choice that we will later make, 
and that seems to be particularly favorable, will violate this condition [ef. Sec- 
tions 7-9, in particular (53)]. 

Instead of working with the coefficients a,, themselves, we can also work with 
the corresponding polynomials 


k 
(14) PZ) = > ay Z' (k = 0, 1, 2,---). 
l=0 


Then (13) becomes 
(15) P,(1) = 1. 


Thus P;(Z) is a kth order polynomial fulfilling (15), and (so far) subject to no 
other restrictions. 
Now (12) becomes, using (10) [and (15)], 


(16) {n‘, a} = P.(E) {&, a}, 
or equivalently, using (11), 
(17) n‘ — &* = P,(G) (E — &). 


The relationship between (16), (17) is similar to that between (10), (11), as dis- 
cussed immediately after (11). 

The broader convergence problem for the iterative-and-mean procedure is this: 
Choose the a,:, i.e., the sequence [Po(Z), Pi(Z), P2(Z),---], [Px(Z) a kth order 
polynomial fulfilling (15), ef. above], so that 
(18) lim n‘ = & 


ko 
for all choices of the starting point &. (18) can be also be written like this: 
(19) lim D(n* — &) = 0, 
kw 


where D(€) is any norm in the space of all mth order vectors & (i.e., in n-dimen- 
sional Euclidean space), which is equivalent to the ordinary (Euclidean) topology 
of that space. We will make a specific choice of D(€) soon: Section 7 (30). 

The ordinary iteration convergence problem (without means, cf. the beginning 
of 4) corresponds to the choice P,(Z) = Z* (k = 0, 1, 2,---) for the sequence 
[Po(Z), Pi(Z), Po(Z),-- +}. 

5. We will now consider the broad convergence problem (iterative-and-mean 
procedure, cf. 4 above) in more specific detail. 
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(19) works with 


(20) di. - D(n* -= &*) (k = 0, * 2,--- ), 

and its requirement is 

(21) lim d, = 0 (for all &). 
k>oo 


Hence the relevant quantity is d, , and we must concentrate on estimating its size 
(for all &). 
Combining (20) and (17) gives 


(22) d, = D[P;i(G)e, 
where 
o =f — &. 


(22), (21) show that the convergence problem is actually one of the convergence 
of the matrices P;,(G)(k — «) to zero. 

Thus the problem presents itself in this form: Given a matrix G, what condi- 
tions must a sequence of polynomials [Po(Z), P:(Z), P2(Z),---] fulfill so as to 
have 


(23) lim P;.(G) = 0. 
k>w 


The answer is well known: Let A; (7 = 1,---, uw, of course » < n) be the char- 
acteristic roots of G, and let e; (= 1, 2,---) be the order of the elementary divisor 
of G that corresponds to ;. Denote the pth derivative of P.(Z) by P,”’(Z). 
Then the necessary and sufficient condition for the validity of (23) is this: 


(24) lim P, (A;) = 0 
k>w 


for all those combinations 7(= 1,---, u), o (= 0,1,---) for which e; > p. When 
all e; = 1, then (24) becomes simply 


(25) lim P;.(A;) = 0 (¢ = 1,--+, mw). 
ko 


For a Hermitian G, in particular, this is always the case. 

We saw in 4, that the P,(Z) are subject to the condition (15): P,(1) = 1. 
Hence (24), i.e., (23), is unfulfillable if some A; = 1, i.e., if 1 is a characteristic 
root of G. In other words: The condition (8c), i.e., the strong condition of 3, is 
reimposed for this reason. [This could have been seen directly too: If that condi- 
tion fails, then for some & ¥ &* there is &' = &, hence all &* = £, hence all n‘ = &, 
and so limy.«. n' = & ¥ &*, contradicting (18).] 

If, on the other hand, (8c), i.e., the strong condition in 3, holds, i.e., if all \; + 
1, then it is not difficult to see that (24) and (15) are compatible. Indeed, even a 
fixed P(Z) {for all P.(Z) with k = the precise order of P(Z), which is }-; e’, ef. 
below] will do: 


P(Z) = el ](Z — i)", 
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with c determined from (15), meets all requirements. This, however, is of small 
practical importance, since the A; may not be known, and the above expression for 
P(Z) may in any case be too complicated for actual evaluation. If it is only known 
that the A, lie in the interior of a certain (bounded and closed) domain A, then a 
sequence [Po(Z), Pi(Z), P2(Z),---] of the desired kind can still be specified, if 
(and only if) A does not separate 1 from ©. We will, however, not go here into this 
matter any further. 

6. The ordinary iterative procedure corresponds, as we observed at the end of 
4, to the choice P,(Z) = Z*. Hence P,”(Z) = k(k — 1) --- (k — p + 1)Z"”. 
Therefore the convergence criterion (24) requires precisely, that all | A;| < 1. We 
state this explicitly: 


The ordinary iterative procedure converges (cf. the beginning of 4) | 
(26) ; 
if and only if | | < 1 for all characteristic values \ of G. } 
As we saw in the last part of 5, this condition is by no means necessary for the 
convergence of some suitable iterative-and-mean ‘procedure. We will nevertheless 
limit ourselves to this case: 


(27) |X| <1 for all characteristic roots \ of G. 


In addition, we will assume that G is Hermitian, because this covers certain 
important applications, and permits the employment of some rather effective 
methods. We restate this: 


(28) G is Hermitian. 


(28) implies that all characteristic roots (or characteristic values) of G are 
real. Hence (27) becomes this: 


(29) —1<X <1 for all characteristic values X of G. 


Under these conditions the ordinary iterative procedure, i.e., the choice (20), 
P,(Z) = Z* (ef. above), is adequate, i.e., it guarantees convergence. We wish, 
however, to determine that iterative-and-mean procedure, i.e., that sequence 
[Po(Z), P:1(Z), P2(Z),---] for which this convergence is (uniformly) fastest. We 
will, therefore reconsider the convergence problem under this aspect, subject to 
the restrictions (28), (29). 

7. We want to choose the sequence [Po(Z), Pi1(Z), P2(Z),---] so as to obtain 
the uniformly fastest possible convergence. This convergence is to be taken in the 
sense of (21), i.e., we want to make for each k the d, of (20) as small as possible. 
{By (22) this means that we want to make D[P;(G@)«] as small as possible.} This 
should be true, in some suitable sense, uniformly—i.e., uniformly in the variables of 
(22). These variables are [since k is given and P;(Z) is being looked for] G and 
w(= & — &*). Let us therefore examine the meaning of uniformity with respect to 
G and w. 

First, since we are now dealing with a situation in which a Hermitian matrix, 
G, occupies a central role, it is reasonable to prescribe that the norm D(£&) be the 
Euclidean norm 


(30) D(§)= Vie ’ 
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(& ,--+, & are the [complex, numerical] components of the [nth order] vector é). 
(Cf. the remark following (19) in 4) 

It is also useful to introduce, for a general (nth order) matrix K, the concepts 
of the ‘‘upper bound”’ | K |, and of the “lower bound” | K |; : 








(31a) | K je = Max 2¢%®) _ win (C such that D(KE) < CD(®)I, 
exo D() 

(31b) | K |, = Min 2°X®) _ ax [C such that D(KE) = CD(®)). 
exo D() 


Second, let us consider the variability of G._G is subject to the requirements 
(28), (29), i.e., it must be Hermitian and fulfill (29). It is clearly logical to make 
the requirement (29) with uniformity, i.e., to require for a suitable e > 0 that 


(32a) -(l-e)SASl1l—eE 
for all characteristic values \ of G, or equivalently 
(32b) lIG|..sl—e. 


Third, let us consider the variability of . We want to make D[P;(G@)o] as 
small as possible (cf. above). w» = 0 is uninteresting, and it is plausible that we 
should want to make the ratio D[P;(G)@]/D(w) uniformly small for all o, i.e., to 
make 


Max PUPs(@)e 
* #0 Do) 


small. This means, by (3la), that we want to make | P;.(G@) |, small. 

Combining the second and the third remarks, we see that we want to make 
| Pi(@) |, uniformly small for all Hermitian G that fulfill (32a) [ie., (32b)]. That 
is, we want to minimize 
(33) a = Max (| P:(@) |.)= Min{C such that for all w and all Hermitian 

7 “\elesi-e G with |@|. < 1— «, D[Pi(G)o] < CD(o)}. 

Now | Pi.(G@) |. is the maximum P;(A), where runs over all characteristic 
values of G. In view of the equivalence of (32a) and (32b), the precise limitation 
on these \ is —(1 — e) SA S 1 — ec. Hence the first part of (33) can be rewritten 
(34) a = Max (| Pe(A) |). 

—(l—e) SAS (1—e) 
Thus we are looking for that kth order polynomial P;,(Z), fulfilling (15), 
P,, (1) = 1, for which 
Max = (| Px(Z) |) 
—(l—e) S$ Zs 1-—e€ 
is minimal. Equivalently: We are looking for that kth order polynomial Q,(Z), 
fulfilling | Q.(Z) | S 1 for all Z in —(1 — €) S Z S (1 — e€), for which Q,(1) 
is maximal. Indeed, for this Q;(1) clearly 


Max = (| Q(Z) |) = 1, 


—(l-—e)<Z<1 
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and 

3 _ Q(Z) 

(35) PAZ) = Oo 

" = 2 << 
(36) «ae P,(Z)|) = Q.(1) ' 


Again equivalently: We are looking for that kth order polynomial R,(Z), fulfilling 
| R(Z) | S 1 forall Zin —1 S Z S 1, for which R,{1/(1 — «)] is maximal. Clearly 


Z 
R.(;4 .). 


8. The last problem in 7 [the one relative to R,(Z)] is classical. It has been 
solved by Chebyshev [2]. The R,(Z) in question is the kth Chebyshev-polynomial, 
defined by 


(38) R, (cos u) = cos (ku). 


(37) Q:(Z) 


It is clear from (38) that R.(Z) is the kth order polynomial, and that —1 < 
Z = 1 implies | R.i(Z) | = 1, as desired. Putting u = i gives R,i(Ch v) = Ch(kv), 
putting ¢” = x gives Ryl}(z + 2”)] = 3(2* + 2“), and putting = Z+- V2 —1 
gives 


(39) RZ) = 3[((Z2 + JZ -— 1) + (24+ V7 —-1)"I. 


Now putting Z = 1/(1 — e) gives 


. 7, we a —k 
(40) r( y= [44S ey +[i+ve ] \. 
—e 2 1—e l—e 


Combining (37) with (35), (36) gives: 


R.( . ) 
(41) P,(Z) = A—< 




















Ms PES RES Ot enngcnes.. 
(42) anna (\P.(Z)I) ( 1 ) 
R. | —— 


It is worthwhile to compare the efficiency of this scheme with that of the ordi- 
nary iterative procedure (without means), i.e., with the choice P,(Z) = Z* (cf. the 
end of 4 and the beginning of 6). 

Consider first the present choice for P;.(Z) [i.e., (41)]. The logarithm of the 
first term in the bracket on the right hand side of (40) is 


In ( + f(2 — i) ot 


1—e 





i.e., for ¢ K 1 it is ~ ~/2e-k. The logarithm of the second term is correspondingly 
~ —4+/2e-k. Assume furthermore +/2e-k > 1, then the first term is dominant, 
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' 
1.€., 
| 
Ry (, —_ .) ~ 3c" (46 
with hi ~ +/2e-k, and so by (42) | wot 
(43a) Max (| Pi(Z)|) = 20" 
—(1—e) $ ZS 1-e 
with hy ~ +/2e-k, if ¢ K 1, VW2e-k > 1. | (47 
Consider next the choice P,;(Z) = Z*. Then clearly 
Max (| P.(Z)|) = (1 — e)*. Put 
—(l—e) < Z<1-« 
The logarithm of the right hand side is In(1 — e)-k, ie., for «eK 1 it is ~ e-k. (48 
Hence in this case 
(43b) Max (|P.(Z)|) =@€”" 
—(l—e) < Z<1-—e i (49 


with he ~ e-k, if e <1. 

Comparing (43a) and (48b), and remembering (34) shows that the speed of 
uniform convergence, i.e., the speed of decrease of d, , compares as follows for the (50 
choices of P;.(Z) under consideration—namely, the “optimum” choice of P;(Z) 
fie., (41)], and the “ordinary” (no means!) choice of P.(Z) (i.e., = Z*): In the 


first case the increase of k that e '-folds d, (asymptotically!) is Ak = 1/+/2e, in 











the second case that increase is Ak = 1/¢. Thus the first choice accelerates the con- (51 
vergence over the second choice in the ratio ~/2e:« = +/2/e. 

9. Let us now return to the definition (38) of R.(Z) [on which the “optimum” 
definition (41) of Pi(Z) is based]. (38) is transcendental, the equivalent (39) is (59 
irrational. It is desirable to replace these by a rational definition. Such a definition 
obtains, in the form of a two-step recursion, from the identity (52 

cos [(k + 1)u] + cos [(k — 1)u] = 2 cos u cos (ku). Ne 
In view of (38) this gives 

Rigs(Z) + Res(Z) = 2ZR,(Z), 

1.€., Ow 
(44) Reai(Z) = 2ZRi(Z) — Ria(Z) (k = 1, 2,---). 
This relation, together with the “starting conditions”’ her 
(45) R(Z)=1, RZ) =Z, (53 
defines the Ri.(Z) completely. Fis 

Now (41) permits us to pass to P;,(Z). Then (44) becomes 

1 1 
Z Ra (; = :) ™ (; = :) 
Prii(Z) = 2 / P,(Z) SF eee P,-1(Z), her 








l—e 1 1 
Ress (; —- :) Ress (; _ :) 








ws 
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~I 
or 


1.e., 





(46) Piss(Z) = 2Z — PAZ) — diss: ax Prs(Z), 
- = 


where 

Ra(; 1.) 
(47) er a 
7 (; - :) 


Putting Z = 1/(1 — e) in (44) and dividing by R,{1/(1 — e€)] gives 











(48) =p. 
Through (41), (45) becomes 

(49) PZ) = 1, P(Z) =Z 
Also, (45) gives 

(50) a4 =1-e«. 
It is convenient to introduce 

(51) b= 
Then (50), (48) give 

(52a) b = 1, 

(ms 1 

(52b) bw = 5-1 - oh, (kK = 1,2, ---). 


Next, (46) gives 
Pryi(Z) = 2igsZ Pi(Z) — (1 — €)” besa de Pra(Z). 
Owing to (52b) 
(1 — ©)? bear be = Zeya — 1, 
hence the above equation can also be written like this: 
(53) Prai(Z) = 2eys (ZPi(Z) — Pra(Z)] + Pra(Z) (k = 1, 2,---). 
Finally by (34), (42) 





hence by (45), (47) 
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and so by (51) 
(54) d. = (1 — €)* bh, --- &. 


10. We can now pass from the P;(Z) to the n“, of course with the help of (16). 
We replace Z by the 2nth order matrix EF in both equations of (49) as well as in 
(53). Thus in all three equations both sides become 2nth order matrices. We apply 
these to the 2nth order vector {&, a}. In this way three equations obtain, each one 
having 2nth order vectors on both sides. These are as follows: 

From the first equation of (49), using (16): 

{n', a} = {&, al, 
1.€., 

From the second equation of (49), using (16): 


1 


{n, a} = EE, a}, 
i.e., recalling (2), (3): 
(56) n' = F{E, a}. 
From (53), using (16): 
{nk a} = 2bear(Ef{n', a} — {ns a}) + {n, a}, 

i.e., again recalling (2), (3): 
{n‘* a} = Qbisal( Fin‘, a}, @) — {né, a}] + {n‘", a} 

= esi Fin', a} — n") + {né”, al, 
1.e., 
(57) ae! = Qbisa(F{n', a} — wh *) + ah, (k= 1,2,---). 


11. We have obtained an inductive definition of the sequence n’, n’, n’, - - - 
This is based on another, inductively defined, (numerical) sequence b, , bo, - - - . 
Actually the two inductions can proceed concurrently. We will now restate these. 

The 6, induction is given by (52a), (52b): 


(Ia) b, = 1, 

(Ib) — sae, (k = 1,2,- ++) 
The n‘ induction is given by (55), (56), (57): 

(IIa) n’ = £, 

(IIb) n' = Fi, a}, 

(IIc) no? = Qbei(F{n*, a} — nS") + (k= 1,2,---). 


We also restate the formula (54) for d : 


(III) dk = | ™ e)‘by "a? by . 
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This can be expressed inductively: 
(IIIa) dy = 1, 
(IIIa) , = (1 —= € beards (k = 0, - 2, of #'S 33 


It is worthwhile to compare this process, and in particular its central piece (I1) 
(which produces the sequence n,n,n,-- -), with the ordinary iterative process, 
i.e., with (9) (which produces the sequence &’, &', &’, - - -). (II) and (9) give the 
same n‘ and z* for k = 0, 1, but they differ for k = 2, 3,- - - , ie., for n** and 
r‘*" for k = 1, 2, - - - . Even here the first step in forming n‘** is the same as the 
(only) step in forming é*~’, the application of the original correction step F{- - - , a} 
(ef. 2). In forming n‘*’, however, this is followed by the further step 2bisi (- - - 
— vn") + a’. This is clearly an extrapolation from n’' with the (excess) factor 
2bi41 — 1. Note, that (I) implies 3 < b; < 1 (foralll = 1,2, ---), henceO < 2b, — 
1 < 1. Thus the extrapolation (excess) factor lies between 0 and 1. 

Now it is by no means unusual that an iteratiye correction method is improved 
by combination with extrapolation steps. The noteworthy circumstance is rather, 
that, in going from n‘ to n‘*’, the extrapolation should issue from n‘~*. It is also of 
interest, that a “universal” and “optimum” sequence of extrapolation factors 
(i.e., the 2b.4; — 1) could be determined [by (1)]. 

12. The procedure summarized in 11 is complete, but it is based on the knowl- 
edge of a Hermitian G fulfilling (32a) [or equivalently (32b)]. Thus there remains 
the problem of constructing such a G. 

More precisely, we need the F of (2), i.e., the G, H of (4). These are linked by 
the relation (7), which we restate: 


(58) G = 1 — HA. 
A is, of course, given. Thus H is arbitrary, it determines G by (58), and this G must 
then be Hermitian and fulfilling (32a). These conditions can also be stated in terms 


of 1 — G = HA: The Hermitian character of G is equivalent to that of HA. (32a) 
is equivalent to 


(59) eSAS2-¢e 
for all characteristic values of \ of HA. 
We repeat: We are looking for an H that makes HA Hermitian and fulfills (59). 
We will now describe two procedures that achieve this: 
First, put 
(60) H = aA* (a > 0).T 
Then (58) gives 
(61) G =I — aA*A. 
HA = aA*A is obviously Hermitian, it is also positive-definite. Hence the smallest 


characteristic value of HA is | aA*A |, = a| A*A |, = a(| A |:)*, and the largest 
characteristic value of HA is | aA*A |, = a| A*A |, = a(| A |.)°. Consequently 











+ A* is the “‘adjoint”’ of A, i.e., its complex-conjugate transposed. 
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(59) means that 


(62a) a(| A |:)* 2 «, 
(62b) a(|A|.)><2-—e 

Assume that we know that 
(63) 0<aZ|A|,5|A|l, SB, 

i.e., so that a, b are known. Then (62a), (62b) can be guaranteed by prescribing 
(64) aa=«¢ ab=2-<., 

i.e., 

(65) a= sy 

(66) e= a 

Now put 
(67) fm?. 

Then (66) becomes 
2 
(68) e= P +i . 

Second, assume that A is Hermitian and positive-definite. In this case put 
(69) H =al (a > 0). 
Then (58) gives 
(70) G =I — aA. 


HA = @A is clearly Hermitian and positive-definite. Hence the smallest charac- 
teristic value of HA is | aA |; = a| A |,, and the largest characteristic value of 
HA is | aA |, = a| A |, . Consequently (59) means that 





(71a) a| A | Po €, 


(71b) ajlAl,s2—e. 


Assuming again the validity of (63), we can guarantee (71a), (71b) by prescribing 





(72) aad = €, ab = 2-—e, 
1.€., 
2 
(73) lies a+b’ 
(74) = 2a 








(7 


0). 


» of 


ing 
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Using (67) again, (74) becomes 
—— 

f+1’ 

13. The results of 12 deserve restatement and some comments. The common 
assumptions of both parts of 12 are (63), (67): 


(75) € 





(IVa) 0<a5|A|,5 |All. SB, 
(IVb) re 
a 

The result of the first part is contained in (60), (61), (65), (66), (68): 
(Va) H = aA* 
(Vb) G = I — aA*A, 

. 2 
(\V Cc) a= e+e ; 
(Vd) a aw 


A being otherwise unrestricted. 
The result of the second part is contained in (69), (70), (73), (74), (75): 


(Via) H = al, 
(VIb) G = 1 — aA, 
. oul 
(VIe) a= s+ b? 
, 2a Sti 2 
(¥aa) ets Stel = T 


A being assumed Hermitian and positive-definite. 

(Va), (Vb) show that the first case is related to the iterative “‘steepest descent” 
methods; (VIa), (VIb) show that the second case is related to the iterative “‘re- 
laxation” methods. 

Our derivation makes it plausible why the former are of universal applicability, 
while the latter are limited to Hermitian and positive-definite matrices—.e., if the 
problem arises from the difference equation treatment of partial differential equa- 
tions of the elliptic type [2nd order, s (= 2, 3, ---) variables, cf. 1 and again 14}, 
to the self-adjoint, elliptic case. 

In general f >> 1. Then in the first case « ~ 2f~* [by (Vd)], and in the second case 
e~ 2f* [by (VId)]. Thus the first case gives a much smaller « than the second case, 
i.e., in view of the remarks at the end of 8, a much slower convergence of the itera- 
tive process. 

This observation illustrates the general experience that whenever relaxation- 
type procedures are applicable, the convergence is significantly faster than other- 
wise. 
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14. We now pass to the consideration of the difference equation system for an 
elliptic partial differential equation [2nd order, s (= 2, 3, ---) variables]. 
Let the partial differential equation be 


(76) ~_ ~< (a' 2) =a, 


i=] Ox; 


where 2; , --- , x, are the independent variables, § = £(x: , --- , 2.) is the dependent 
variable, and a’ = a'(z,,---, 2), (¢ = 1,---,8) anda = a(xz,,---, 2,) are 
known functions of x, , --- , x, . Also 
(77) 0<a@ sa'(m,---,2,) <b’ (t =1,---,s8), 
the a’, b' (¢ = 1, --- , s) being known constants. Finally the domain of the 2; , 
x, 18 
(78) Os2,s L; (i = 1, , 8) 
and the boundary condition is 
(79) €=0 for z;=0 or L; (¢ = 1, 8) 
In order to pass to difference equations, we introduce a lattice 
L; 
(80) Xi= ni Az; | Ar; = N. ’ 
where in some cases 
(80a) n =0,1---,N;-—1,N:, 
in others - 
(80b) n=1,---,Ni-1, 
and in others again 
(80c ) mn = $,%,°--,Ni — } («= 1, $) 
Of course, VN; = 2, 3, --- , and it expresses the fineness of this lattice in the 
x;-direction. 
We write 
(81) E(x, ee 2s) = Eni---m ’ 
using (80a). These &,,...,, are the unknowns, but since (79) gives 
(82) £....¢,= 0 for », = 0,N, (¢ = 1, ---, 8), 


the unknown character is actually restricted to (80b). 
It is convenient to use with a’ (80b) for the n; with 7 ¥ 7, and (80c) for 7; : 


(83) a’(a,-+**,2%s) = a 
and for a (80b) throughout: 
(84) (21, °° Le) = Ogy---9, 5 


these being known quantities. 
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Now (76) is best stated for (80b). It becomes 
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; 
Oy ---95—----45 (Eps -ong-o-ne a Ea1---95-1---9)) = Ogi---me - 


We can view (85) as the equivalent of (1), with the following provisos: The 
complex m , --- , m., according to (80b), stands for the vector-index in (1). Hence 


the order of the matrix A is 


(86) n= [I (N; - 1). 


The &,,.--,, are therefore the components of the (unknown) vector &, the a,,...,, are 
the components of the (known) vector a. The left-hand side of (85) then defines 


the (known) matrix A. Hence 
(87) A=>> (*:) Ai, 
i=l L; 


where the matrix A; is defined by 


At =f, 


+ i 
(88) Seam = ay --ngth---me (San---nett---ne — Say---05---9e) 


? 
Fg y-+-nj teeny Ena esengeeeng — Eqye-ng—teong)- 


Furthermore, clearly 
(89) A; = BSB; 
[ef. footnote on page 177, where 
BE=&, 
— a ao 


(90) 


In order to apply the results of 13, we now need estimates of | 


accordance with (IV) in 13. From (87) 


lal (Ni 
|A|:2 pam FL 


i=l 


IV 





) 
) | Ag |2 | 


Al. < > (*) 
| = |e fant a 444i g ‘| 


(91) 





From (89) 
| As |: = (| Bi |2)’, ) 


(92) . 
A; a _ (| B; wf 





Enu---mis-oms ) = 


uy 


in 


Finally, designate A; , B; with aj,....,...., = 1 fef. the remark preceding (83) !] by 


A,’, B?. Then clearly 


& 
IV 


| > Va‘ | B? |:, ) 
(93) 


Voi | Be \u. J 
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Combining both sides of (93) with (92) gives 


| Asli ‘TAS |e, ) 
(94) 


| A; i. s b’ A? ies 


IV 
Q 


and combining (94) with (91) gives 


| - at N; ? " 
[Al z Da‘ (™) [Ae |, | 
i=l i | 
(95) > v2 r 
Ales D8 (M) | as. 
i=l 4G 


Now consider A,’. Applying (88) with a},...,;-..., = 1 shows, that the role of the 
. > 2 . ° - j 0 0; . ° 
n;,j ~ 7, is now irrelevant in determining | A; |;, | Ai |, . Hence we may write in 
place of (88) 
A PE 3 = r 


(96) , , 
En; _ —b541 <5 2&,; _ fy;—1 -) 


| 
_ 


This operator is Hermitian, its characteristic vectors are the —™*(m; = 
N; — 1) with 


- TM; Ni 








(97) gS) = sin v 

the characteristic value of —”* being 

(98) AN = 2 — 2 cos T™ = 4 sin’? T. 
Hence 


|A°|, = Mind™ = = 4 sin’ : 

(99) m; ava 
“gm } si a 

| Ai |u Bios 2 p 4 cos ON.'| 


bo 





Combining (95) and (99) gives 
1 : de N; P 202 T 
|A|,2 4 > a; (*) sin’ | 
= = N; , 2 fF [ 
| A le 4 Y b; (*) cos ON P 


Hence we can put in (IVa) 
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Therefore (IVb) gives 











(102) == 5 
a N; 2 
% a{ ¥) sin —— 
i=l L; QN 
If 
(103) Max = = M Max N; = N, 
i=l,---,s G; i=l,---,8 
Then (102) gives 
104) < M ot? me . 
( f=Me oN 
and, since 
k ap © 2N 
‘non =on’ “ne, 
a fortiori 
‘ 4— 
T 


Now (VId) in 13 gives « = 2/(f + 1) and the remarks at the end of 8 give for 
the error-e '-folding increase of k (asymptotically!) 





wm ten ypet~ iyi, 
i.e., 
(106) ak ~ 5 Vj. 
Hence, in view of (105), 
(107) ak <* Vin. 


15. We restate the results of 14. The elliptic partial differential equation is given 
in (76), the subsidiary conditions in (77)—(79), the lattice is defined in (80) and 
(81)—(84), the difference equation system in (85). We do not restate these. 

The a, b, f of (IV) in 13 are given in (101), (102): 








r = “ N; , - 2 F&F 
(VIIa) a=4>a(™) sin ON” 
. _,oz (NY 2 & 
(VIIb) b= 45 (™) 00s" Sy. 

= 2 N; 2 fF 
5 ( ) cos 
(Vile) ja ee 
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If 
(VIIIa) 
then 


(VIIIb) 
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The Determination of the Chebyshev 
Approximating Polynomial for a 
Differentiable Function 


F. D. Murnaghan and J. W. Wrench, Jr. 


If f(z) is continuous over any interval, which we may take, without loss of 
generality, to be the interval —1 < z < 1, there exists a unique polynomial P,, *(x), 
of given maximum degree n, which is such that the maximum of | f(z) — P,*(zx) | 
over —1 S x S 1 is less than the maximum of | f(x) — P,(x) | over —1 S x <1, 
where P,(2) is any other polynomial of degree not exceeding n. The polynomial 
P,*(x) is known as the Chebyshev approximation, of maximum degree n, to f(z) 
over —1 S 2 S 1. It is characterized by the fact that f(z) — P,*(x) assumes 
extreme values at n + 2 points, at least, of the interval —1 < x < 1, these ex- 
treme values being equal in magnitude and alternating in sign [1]. We refer to the 
points of any such set of n + 2 points as critical points, and we denote them by 


* * . 
(ay*, --- , Un42), Where —1 S a,* < --- < Xny2 S 1. Thus, the end points, +1, 
of the interval —1 < x S 1 may be critical points, but at least n of the n + 2 
critical points, namely, 22*, --- , 2x41, are interior points of this interval. 


We assume that f(z) is not only continuous, but also differentiable, over —1 < 
x < 1, and so the derivative of f(x) — P,*(x) is zero at each of the n points x.*, 

- ,tng1 . If the derivative of f(x) — P,*(x) cannot be zero more than n times, it 
follows that x,* = —1, 2.42 = 1 and that the derivative of f(x) — P,*(x) has pre- 
cisely n zeros that are interior points of the interval —1 S zx < 1. 

The polynomial P,,*(x) is an odd function of x when f(z) is odd, and is an even 
function of x when f(z) is even. Thus, when f(z) is odd we may take n to be even, 
the maximum degree of P,,*(x) being n — 1, and when f(z) is even we may take n 
to be odd, the maximum degree of P,,*(x) being again n — 1. In these cases the 
critical points are distributed symmetrically about the mid-point z = 0 of the 
interval —1 < x < 1, and we may confine our attention to the part 0 < x < 1 
of this interval. When f(x) is odd, so that n is even, the number of critical points is 
even and xz = 0 is not a critical point; on the other hand, when f(x) is even, the 
number of critical points is odd and x = 0 is a critical point. When f(x) is odd, or 
even, and x = 1 isa critical point, we change our notation and denote the positive 
interior critical points by 1* < 22* < --- < 2,*, where n = 2k in the first case, 
and n = 2k + 1 in the second. For example, when f(x) = are tan 2, P(x) is an 
odd polynomial of degree £2k — 1, and so the derivative of are tan x — P(x) 
cannot vanish more than 2k times; this implies that the points +1 are critical 
points, and, in addition, since this derivative must vanish 2k times, that P(x) is 
of degree 2k — 1. Similarly, when f(z) = cos mz, m > 0, Pr41 is an even poly- 
nomial of degree 2k, the points 0 and 1 being critical points. 

It is clear that P,*(2x) is easily determined if any set 1* < m* < --- < Te42 
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of critical points is known; indeed, the n + 2 equations f(x,.*) — Pna*(a*) = 
(—1)*"w, k = 1, --- ,n + 2, constitute a set of n + 2 linear equations for w and 
the n + 1 coefficients of P,,*(x). If 71 < a2 < --+ < Xn42is any set of n + 2 points 
of the interval —1 < xz < 1, and we write f(a.) — P,(a) = (—1)*"E, k = 1, 

- ,n + 2, these equations determine EF and the n + 1 coefficients of P,(x), and 
the function EF of the n + 2 variables (2, --- , tn42) has an absolute maximum 
at (a:*, --- , 2242). Thus, the derivative of E with respect to each of the n variables 
2, °°* , tnqi is zero at (m*, ---, Sese)s and this implies that the derivative of 
each of the n + 1 coefficients of P,(2) with respect to each of the n variables 
Ze, °** » Tn41 18 Zero at (7*, ---, Ta+2). Hence, these coefficients are insensitive to 
small changes of the variables (12, --- , 2,41) when these variables have the 
values 22*, ---, tey1 and 2 = 2%, Inge = La4e- 

The method by which we determine P,,*(x) is an iterative one. Let us suppose 
that the points +1 are critical points, so that there are n interior critical points, 
which we denote, changing slightly our previous notation, by 1* < x* < --- < 
z,*. Let us suppose further, that we are in possession of a polynomial, P,° (x), of 
degree <n, which we term our entering polynomial and which possesses the follow- 
ing property: The difference f(x) — P,,° (x) assumes extreme values of alternating 
sign at n + 2 points of the interval —1 < x S 1. We denote by a, °** she 
approximations to the second, third, --- , (n + 1)st of these points and we regard 
a, +--+, 2° as approximations, in the first cycle of an iterative procedure, to 
ai*, «++ , tn*. We deterniine the approximating polynomial of degree <n, P,,” (2), 
with which we end the first, and begin the second, cycle of this procedure by means 
of the n + 1 linear equations obtained by eliminating E™ from the n + 2 linear 
equations 


f(—1) —_ P,(-1) - ge. f(a) tad P,” (x) - (—1)°2", 
k=1,---,n;  f(1) — P,(1) = (—1)"°"E". 


Denoting f(—1) — P,(—1) by oo, f(a”) — Pa (ae) by &”, k = 1, +--+, 0, 
and f(1) — P,(1) by 5.7, we can write these n + 2 equations as 


6P,(—1) a a.” * E”. bP, (ax?) ah = (—1)"E", k = 1, ae nN; 
bP. (1) = bny41 — (—1)"7B® 


where 5P,,° (x) denotes the polynomial, of degree <n, P,? (x) — P, (x). E™ is 
conveniently eliminated by combining the last n + 1 of these n + 2 equations al- 
ternately by addition and subtraction with the first, and the n + 1 coefficients of 
éP,,(x) are obtained by solving the resulting n +1 linear equations. Then the 
coefficients of P,,"’(x) are obtained by adding each of the coefficients of 5P, (x) 
to the corresponding coefficients of P,,° (x). 

The first step in the second cycle of the iterative procedure is the determination 
of new approximations x,”, --- , x,” to a*, «++ , an*. Just as 4", --- , an” were 
approximations to the zeros of D[f(x) — P,’(x)], where D denotes differentiation 
with respect to x, so 2”, --- , 2» are approximations to the zeros of D[f(x) — 
P, (x)]. Writing x = a + da, k = 1, «++, , we see that the value of 
Dif(x) — Px (x)] at ae + 62, must be the same as the value of D[6P,(x)] at 
a + da,"”, and this is the same, to the first order of infinitesimals, as the value 
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of D[6P,, (x)] at zx. Thus, to the first order of infinitesimals, the value of Dif(x) — 
P,”(a)) at 2 plus the value of D*[f(z) — P,(x)] at 2,” times dx,” is equal 
to the value of D[éP,°(x)] at 2°, so that 5x," is the negative of the quotient of 
the value of D[f(z) — P,“(x)] at 2,” by the value of D’'{f(z) — P,(x)] at x". 
We denote the value of D[f(xz) — P,"(z)] at x” by a”, k = 1, --- , n, so that 
dx,” is the negative of the quotient of «” by the value of D’[f(x) — P,(z)] at 


ay”. If «” is zero, x” = 2,” and D{f(z) — P(x) is zero at x = x. In 
order to complete the second cycle, we calculate the n + 2 numbers dp)” = f(—1) — 
Pa'(—1), be = f(t) — Pa (ae), (ke = 1, «++, m), Omar = f(1) — Pa (1), 
and determine, as before, the coefficients of 6P,,°’ (x). If, to the number of decimals 


we are using, the n + 2 numbers 4, --- , 6.4; are equal in absolute value and al- 
ternating in sign, the coefficients of 6P,,"’(x) are all zero and the approximating 
polynomial, P,,” (2), with which we end the second cycle is the same as the ap- 
proximating polynomial, P,,"’(x), with which we began it. 

In a previous publication [2] we determined the numbers 2; , --- , 2, in each cycle 
by solving the equation D[f(z) — P,(x)] = 0, where P,(z) is the approximating 
polynomial, of degree <n, with which we begin the cycle, but the less exacting 
method of the present paper is equally effective. 

It remains only to describe the selection of our entering polynomial approxima- 
tion, P, (x), of degree < n, and the determination of the approximations z,"”, - - - , 
zt,” to the n points of the interval —1 < x < 1 at which Dif(z) — P,(z)] is 
zero. On setting x = cos 6, we see that f(x) becomes a function, F(@), of @ defined 
over 0 S 6 S z, and we write the Fourier cosine series of F(@) as 4a) + a, cos @ 
+ a, cos 26 + --- . Then cos mé@ is a polynomial function, 7,,,(z), of x of degree 
m, which is known as the mth Chebyshev polynomial, m = 0, 1, 2, --- , and 
da) + a7Ti(x) + a2T2(x) + --- is known as the Chebyshev expansion of f(z). 
The sum of the first n + 1 terms of this Chebyshev expansion of f(x) is a poly- 
nomial function, of degree <n, of z, and it is this polynomial function that we take 
as P,, (x). We say that P,, (x) is furnished by the truncated Chebyshev expansion 
(the truncation taking place at the term which involves 7,,(z)). Now f(x) — 
Pa (x) = Qn41T n4i(2) + --- , and we take as our approximations to the n points 
of the interval —1 < x < 1 at which Dif(z) — P,,(x)] = 0 the n points of this 
interval at which D[T,4:(x)] = 0, it being assumed that a,,, ~ 0. (If f(x) is odd 
its Chebyshev expansion is of the form a7;(z) + a;7;(z) + --- andn = 2m is 
even; then we truncate this Chebyshev expansion at the term involving T2,_:(2), 
and we take as our approximations to the m points of the interval 0 < x < 1 at 
which D[f(x) — P,, (z)] is zero the m points of this interval at which D[T2m4:(x)] 
is zero. Similar remarks apply to the case where f(x) is even and n = 2m + 1 is 


odd.) Since 


sin (n + 1) 6 


D [Tn41(2x)] = (n+ 1) a6 


, 


we have 


ni” = cos (x - **), k=i,--+ a. 


Example 1. f(x) = are tan z,n = 6. 
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There are three positive interior critical points, which we denote by 2;*, 72*, 23*, 
and 
Pe (x) = 0.9949493662 — 0.2870606362*° + 0.0789371762°, 
since the Chebyshev expansiun of arc tan z is 
3 5 
2 | prix) = 5 T3(2) + e T;(z) “mm cc ‘|, 
where p = 2} — 1 = 0.414213562, to 9 decimals. The first-cycle approximations to 
2*, r2*, 23* are 
= 0.222520934, x2” = 0.623489802, x;"” = 0.900968868, 
and the polynomial approximation with which we end the first cycle is 


Ps (x) = 0.995383022x — 0.2887004402* + 0.0793133072°, 


the values of are tan x — P,, (x) at the points 1,", x2”, 2;"", 1 being 
6, = 0.000676851,  &° = —0.000604555, 
5; = 0.000546760, 4," = —0.000527744, 


respectively. The corresponding results for the second cycle are 
» = 0.205422893, 2” = 0.593832571,  2;° = 0.888813502, 
= 0.000603543, &” = —0.000619441, 
6; = 0.000607728, 4° = —0.000597725, 
and ‘ 
P.” (x) = 0.995357994x2 — 0.2886904172* + 0.0793391732°. 
In the third cycle we find 
= 0.205218790, 2° = 0.593469973, x; = 0.888196372, 
5° = 0.000608588, 6° = —0.000608590, 


6; = 0.000608612, 6° = —0.000608588, 


and 
P. (x) = 0.995357955« — 0.2886902382° + 0.0793390412°5. 
Finally, in the fourth cycle we obtain 
= 0.205219373, «x = 0.593470162, _—z;“” = 0.888196289, 
5: = 0.000608595, 6 = —0.000608595, 4,” = 0.000608595, 
6, = —0,000608595, Pes (x) = Ps (x) 
Thus, to seven decimals, 


Pe*(x) = 0.9953580x — 0.2886902zx° + 0.0793390z°, 
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and the maximum of | are tan z — P,*(x) | over —1 S x S 1 is 0.0006086. Hastings 
[3] has given P,*(x) for are tan zx to six decimals as 0.9953542 — 0.2886792° + 
0.0793312°. 


+2 10° + 1 
e- 


Va Gg a oT = —_ = 
Example 2. f(x) = log peng ior = 1’ mn = 4, 


The number a must be greater than 1; we use the value indicated in order to 
check the work of Hastings. Setting = (a + x)/(a — x), the polynomial P,*(x) 
which we determine will be an approximation to log — over the interval 10° < ¢ < 
10°. 

The Chebyshev expansion of log (a + x)/(a — z) is 


3 5 
4M | prix) + y T(x) + E T(z) + |, 


where M = log e = 0.434294482, to 9 decimals, and p = a — (a — 1)! = 
0.280130000. The entering polynomial approximation is P,’ (x) = 0.448447982z + 
0.0509168942°, and our first-cycle approximations to the two positive interior critical 
points are 


a” = 0.309016994; x2” = 0.809016994. 
The values of log (a + x)/(a — x) — P, (x) at the points x;"", x2", 1 are 
6:” = 0.000572828, 62” = —0.000607958, 6," = 0.000635124, 
respectively, and P,"(x) = 0.4483493552 + 0.0510513052°. In the second cycle 
a” = 0.321484228; 2,” = 0.821954759 
6°” = 0.000600487, 6° = —0.000602901, 6; = 0.000599339 
P.? (x) = 0.448347007x + 0.0510517662", 
and, in the third cycle, 
a = 0.321320097; = x2” = 0.821455202 
6, = 0.000601227; 6° = —0.000601233; 4; = 0.000601227 
P(x) = 0.448346999x + 0.0510517712", 


so that, to seven decimals, P,°(x) coincides with P,° (x), and P,*(x) = 
0.4483470x + 0.05105182°, the maximum value of | log (a + x)/(a — x) — P,*(z) | 
over —1 S x & 1 being 0.0006012. 

If x is replaced by a(x — 1)/(x + 1), there results the approximation 0.8630458 
(2 — 1)/(a + 1)] + 0.3641410[(2 — 1)/(a + 1)]° to log x over the interval 
'<2x< 10. Hastings [3] gives as the coefficients in this approximation the 
numbers 0.86304 and 0.36415, respectively. 

In a recent paper by Barth [4], P.s*(x) for n (a + x)/(a — x), a = (10° + 1)/ 
(10° — 1), is given, to ten decimals, as 


0.8690286986x + 0.27738331952° + 0.2543282307 2°. 








190 F. D. MURNAGHAN AND J. W. WRENCH, JR. 


The correct formula, to seven decimals, for P.*(x) is 
0.86902852 + 0.27738642° + 0.2543195z°, 
the maximum of | In (a + x)/(a — x) — P¢*(x) | over —1 S x S 1 being 0.0000337. 


Example 3. f(t) = In(1+£&),0 S&S 1,n = 4. 

We denote in this example the independent variable by &, instead of x, since the 
interval, 0 < ~ < 1, is not the standard interval —1 < x S 1. The linear trans- 
formation x = 2 — 1 transforms the interval 0 < & < 1 into the interval —1 < 
x S< 1, and the problem of determining P,*() for In (1 + £) over 0 S ~ S 1 is the 
same as that of determining P,*(x) for In 3(3 + x) over -—1 S x J 1. 

The Chebyshev expansion of In $(3 + 2) is 


4 6 
—2 log. 2p + 2 | pr u(x) = E T(x) + a T3(x) |, 


where p = 2' — 1, and our entering polynomial approximation, of degree 4, is 
P,(€) = 0.000069446 + 0.996261948¢ — 0.466442439% + 0.218665484¢° — 
0.055459314é*. Our first-cycle approximations to the four interior critical points 
are 


&° = 0.095491503,  & = 0.345491503, 

&” = 0.654508497, &” = 0.904508497, 
and the values of In (1 + £) — P,(£) at the points 0, &, &, &°, &'°, 1 are 
5 = —0.000069446, 4, = 0.000066650, 6°” = —0.000060948, 

63° = 0.000055988, &° = —0.000053017, 4; = 0.000052055. 


The polynomial approximation, of the fourth degree, with which we end the first 
cycle is, to eight decimals, 


P, () = 0.000059471 + 0.996558114¢ — 0.467864445¢° 
+ 0.220882267¢° — 0.056547698¢". 
In the second cycle we obtain 
£,° = 0.084407707,  &” = 0.318071278, 
&” = 0.629216597,  &° = 0.895475308; 
50” = —0.000059471, 4” = 0.000060703,  &” = —0.000061938; 
8s” = 0.000061315, 8” = —0.000060043, 4” = 0.000059471; 
P,?(€) = 0.000060712 + 0.996540728¢ — 0.467834593¢" 
+ 0.220891205¢* — 0.056571583¢"; 
and, in the third cycle, 
#:° = 0.085058286, &” = 0.319106141, &° = 0.629174645, 


& = 0.895122761; 
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5 = —0.000060712, 4,” = 0.000060716,  &” = —0.000060716, 
5; = 0.000060712, 4,” = —0.000060713, 4; = 0.000060712: 
P;®(€) = 0.000060714 + 0.996540741¢ — 0.467834762¢ 
+ 0.2208915412° — 0.056571768¢". 
Beginning the fourth cycle, we compute 
&° = 0.085060350,  & = 0.319112305, 
é; = 0.629171981,  &,” = 0.895123131, 


and we find that the corresponding values 6; for i = 0, 1, 2, 3, and 4 are all nu- 
merically equal to 0.0000607141, to within a unit in the tenth decimal place. 
Thus, to seven-decimal accuracy in the coefficients, we have 


P.*(—) = 0.0000607 + 0.9965407E — 0.4678348¢’ + 0.2208915¢° — 0.0565718¢". 


The discrepancy between this approximation and the similar one presented in 
our earlier paper [2] is attributable to the premature termination of the iterative 
procedure in that reference, which stemmed from the erroneous belief that the pre- 
cision of the coefficients of the approximating polynomial was comparable to that 
of the maximum difference between that polynomial and the given function. In this 
example the quantities 5,“ have all become stabilized to 10 decimal places, whereas 
the coefficients of the corresponding approximating polynomial are subject to 
errors of approximately a unit in the eighth decimal place. This behavior of the 
coefficients is due to the relatively small value of the determinant of the system 
of equations used for their evaluation. Calculation of such coefficients to ten-place 
accuracy generally will require double-precision operations. 

Hastings [3] gives as an approximating polynomial of degree 4, whose graph is 
arbitrarily required to pass through the origin, the following 


0.99744422 — 0.47128392* + 0.22566852° — 0.0587527z", 


for which the maximum departure from In (1 + 2) over the interval 0 < z S Lis 
0.0000710, in contrast to the value 0.0000607, attained by the Chebyshev approxi- 
mating polynomial of the same degree. 

Example 4. f(x) = cos (r/4)z, n = 3. 

The Chebyshev expansion of cos (/4)z is Jo(4/4) — 2Jo(4/4)T2(x) + 
2J4(4/4)Ts(x) — +--+, where Jx(x/4), for k = 0, 1, --- , is the value at +/4 of 
the Bessel function of the first kind, of order 2k. There is only one positive interior 
critical point z*, and our first-cycle approximation to this is x“’ = cos 4/4 = 
0.707106781. Our entering polynomial approximation is P;° (x) = 0.998068558 — 
0.2928732892x". The values of cos (x/4)x — P; (x) at the points 0, x’, 1 are 


5)? = 0.001931442, 4° = —0.001921422, &” = 0.001911512, 


and P;"(x) = 0.998078551 — 0.2928932192". The coefficient of x* remains the 
same in all the succeeding cycles, so that, in this example, we obtain the coefficient 
of x” in P;*(x) before we begin the second cycle; this simplification is due to the 
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fact that both the points 0 and 1 are critical points, and this implies that the co- 
efficient of x” in P;*(x) is cos r/4 — 1. In the second cycle we obtain 


= 0.705276652, 6” = 0.001921449, 6,” = —0.001921553, 
5.” = 5)”, P(x) = 0.998078499 — 0.2928932192": 
and, in the third cycle, 
’ = 0.705270859, 6° = 0.001921501, 46,° = —0,.001921500, 
i” = 6°, P(x) = 0.998078499 — 0.2928932192°. 


Thus, to seven decimals, P;*(z) 0.9980785 — 0.2928932z", the maximum of 
|cos (4/4)a — P3*(x) | over —1 S x S 1 being 0.0019215. 

We observe in this example that the entering polynomial approximation, P;“’ (x), 
is so good that the maximum of | cos (4/4)a — P;‘(x) |, over —1 S x < 1, is 
0.0019314, which exceeds the corresponding maximum of | cos (7/4)x — P;*(x) 
by less than 0.52 per cent. 


A Il 


Example 5. f(x) = cos (/2)x,n = 5. 


There are two positive interior critical points, 2;* and 2,*, and our first-cycle 
: : (1) « a jo 92/9 TW ‘ 
approximations are 2". = cos r/3 = } and x2” = cos x/6 = 3°/2. The € hebyshev 


expansion of cos (#/2)x is Jo(a/2) — 2Jo(w/2)T2(x) + --- , and our entering 
polynomial approximation, of degree four, is Ps°(x) = Jo(a/2) — 2Jo(4/2)T2(x) 
+ 2Jo(e/4)Ts(x) = 0.999396554 — 1.2227431532° + 0.223936637x*. The values 
of cos (x/2)a — P;(x) at the points 0, 1", x2", 1 are 
5°” = 0.000603446, 6” = —0.000600024, 

52” = 0.000593320, —;°” = —0.000590037. 


Since 0 and 1 are critical points, the coefficients of the approximating polynomial 
P; (a), with which we end the kth cycle (k = 1, 2, ---) satisfy the relation 


2a + gp + 7“ = 1, and this implies that the coefficients of P;*(x) satisfy the 
relation 2a* + 6* + y* = 1. We find that 


Ps (x) = 0.999403304 — 1.2279688027 + 0.223990272z°*. 
In the second cycle we obtain 
= 0.497202761, 2” = 0.864404535; 
do” = 0.000596695, 4° = —0.000596808, 
5.” = 0,000596805, 6; = —0.000596697; 
Ps” (x) = 0.999403231 — 1.2227967332" + 0.223990272z". 
In the third cycle we obtain 
= 0.497195260, 2x = 0.864395279; 
® = 0,000596769, 6° = —0.000596772 
5.” = 0.000596770, 6; = —0.000596770; 


Ps? (x) = Ps” (x) to 9 decimals. 





10- 


~I 


al 


1e 





CHEBYSHEV APPROXIMATING POLYNOMIAL 193 


Hence, we conclude that 
P;*(x) = 0.9994032 — 1.2227967x" + 0.2239903z", 


the maximum of | cos (#/2)z — P;*(x) | over —1 S x S 1 being 0.0005968. 

The iterative procedure described herein for the determination of the Chebyshev 
approximating polynomial P,*(x), of degree not exceeding n, over the interval 
—1 S x S 1, to a given differentiable function f(z) will converge if the difference 
between f(z) and the initial approximating polynomial P,’(z) assumes extreme 
values at n + 2 points of the interval —1 < zx S 1 and if, furthermore, these ex- 
treme values alternate in sign. A proof of this based on the argument of Novodvor- 
skii and Pinsker [5] has been given, and illustrated by a numerical example, in our 
previous publication [2] on this subject. 
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A Method of Computing Eigenvectors 
and Eigenvalues on an 
Analog Computer 


By Lucien Neustadt* 


1. Introduction. Many papers have been written in recent years describing 
methods for finding the eigenvalues and eigenvectors of an arbitrary matrix. 
Most of these methods apply to digital computation, and little attention has been 
paid to methods applicable to analog computers. One such method, which has been 
used successfully in the past, is described in [2]. The analog technique described in 
the present paper has the advantage of yielding both eigenvectors and eigenvalues 
of a real symmetric, or complex Hermitian matrix, and of not requiring a trial 
and error procedure. This is accomplished at the expense of additional computing 
equipment. 


2. Mathematical Formulation. Consider the real, symmetric, n X n matrix 
A = (a;;) with a;; = a;;. We shall state, without proof, the following properties 
of the matrix A [3], [4]: 


(1) There is an orthonormal set of n real eigenvectors e’, e’, --- , e” of the 
matrix A; i.e., there exist n numbers di, 2, ---, An, (the eigenvalues) with 
Ae’ = de’ andi = 1, 2, --- , n; also 

e'-e’ = 3,7 : 


(2) The eigenvalues are all real. 
Suppose that the eigenvalues are arranged in descending order such that 





X\; 2 Ae 2} --- 2 An; then for an arbitrary real, non-zero vector 
oat 
x=/: 
Tn 
n 
(3) A, = sup “—"_____ = sup = sup Ax-x 
x 2 x xx x*x=1 
Dt 


i=l 


The sup is attained at an eigenvector e’, belonging to d,-e' and \; for i > 1 
can be obtained from the same formula if the following added restriction is made: 


x-e’ = 0 j=1,2,°--,t—1. 
e” and X, can also be computed from formula (3) if the sup is replaced by inf. 


Received March 7, 1958, in revised form, August 26, 1958. 
* Now at Space Technology Laboratories, Los Angeles, California. 
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3. Solution Procedure. The technique is based on the analog computer method 
for solving problems in linear programming [5]. This method may be readily adapted 
to problems in nonlinear programming, i.e. to the problem of finding the extreme 
values of a nonlinear function subject to nonliaear restrictions. The problem of 
finding eigenvalues and eigenvectors is just such a problem, namely that of lo- 
cating the extreme values listed above. 

Let x be the radius vector to the point x = (2, ---, 2.) in n space. Let x, 
22, °** , Xn be functions of time. As described in the mathematical formulation, 
e' may be found as follows: 

Let the point z move on the unit sphere until it reaches a steady state corre- 
sponding to a maximum of the function Ax-x. Then e’ equals the steady state of x. 

The vector e* is found in the same way, with the additional restrictions on x 
that x-e' = x-e’ = --- = x-e*"* = 0. 

Also, e” is found in the same manner as e’, except that the steady state now 
corresponds to a minimum of Ax-x. 

In order to find e’, we write a differential equation for the vector x such that the 
steady state solution of the equation is e’. 

Denote Ax-x by f(x) = f(a, 2, +++, 2n). 

If we set x = grad f, the point x will move in such a way as to increase f. To 
insure a steady state with x-x = 1, modify this equation to read 


(1) x = grad f + ka 
where 

lifx-x<s1 

is —lifx-x>1 


and k is a positive constant, chosen such that 
k > 2 max (| A; |, | An |). 


The point z will move as follows: Assume z(0) = 0. The origin is a point of 
unstable equilibrium. Once z moves slightly away from zero it will continue to 
move in a direction which is given by the vector sum of grad f and kx. Because k 
is large enough, x will move toward the boundary of the unit sphere. Once it breaks 
through, again because of the choice of k, it will immediately re-enter. Having re- 
entered, it will immediately break through again, and so on indefinitely. The point 
zx will then oscillate back and forth across the surface of the sphere with a fre- 
quency depending on the rapidity of the switching arrangement which determines 
the sign of e«. Superimposed on this oscillation is a tangential motion produced by 
the tangential component of grad f. The tangential motion will continue until a 
point is reached where grad f has only a normal component. There grad f(x) = ax 
for some constant a. It is easily verified that in general, grad f(z) = 2Ax, so that 
this point corresponds to an eigenvector. Furthermore, at such a point, f(z), where 
z is restricted to lie on the unit sphere, has an extreme value. If x ~ e’, the solution 
is unstable, since any motion of x in the direction of e will be self-reinforcing. A 


stable steady state solution is reached only when x = e. 











196 LUCIEN NEUSTADT 

The eigenvector e”, belonging to the smallest eigenvalue X, , is the steady state 
of the solution to the equation 
(2) x = —grad f + ke. 


If equation (1) is modified, the steady state of x can be made to correspond to 
the eigenvector e’. Consider the differential equation: 


(3) x = grad f + kex + hie: grad f; ;x 
where (a) ¢ and k are defined as above. 
(b) fi(z) = x-e’, 
[1 if f(x) $0 
(c) «1 allt 


\-1. if fi(~) > 0, and 


(d) k; is a large enough positive constant. It is sufficient that k, > 2k. 

The added term will insure that the steady state of x yield a radius vector x 
orthogonal to e’. Since the maximum condition is also satisfied, this steady state 
corresponds to the eigenvector e’. 

The other eigenvectors can be obtained by modifying equations 2 or 3 by the 
addition of similar terms. 

In principle the eigenvalues can be obtained directly from the eigenvectors. 
Given the eigenvector e’, one may compute Ae’. The ratio of the components of 
Ae’ to the components of e’ equals the eigenvalue \;. Since e’ is obtained only 
approximately, this ratio will in general not be constant. 

An excellent approximation to the eigenvalue, knowing an approximate eigen- 
vector e’, is given by - 


Ae‘-e’ 
e’-e? 





[4]. 


4. The Computer Setup. In order to find the eigenvector corresponding to 
the largest (or smallest) eigenvalue of A the following amount of computing equip- 
ment is necessary : 


nm integrating amplifiers, whose outputs are 2 , t2,-°-- ,2%n. 

nm inverting amplifiers, whose outputs are —2,, —%2,°-*-, —In- 

n’ scale factor potentiometers, corresponding to n potentiometers for each 
component of grad f. 

One “switch” to compute ¢. A high gain amplifier with two diodes and two 


voltage sources can be used for the switch. 


Multiplying equipment to compute 2; 72, --- , Zn, a8 well as ex, ef2, °°: , 
ex, . For this purpose n multiplying servos positioned by 2, --- , %, may 
be used. 


One inverting amplifier to compute —. 

In order to find the eigenvector corresponding to dz (or An-1) the following 
additional equipment is necessary: 

A switch to compute « . This can again be a high gain amplifier. 

n potentiometers to compute fi(z). 

n potentiometers to compute the components of grad fi . 
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One inverter to compute —« . 

For every additional pair of eigenvectors more of the same equipment must be 
used. 

The eigenvalues can also be found with the computer. Having computed the 
approximate unit eigenvector e’, \; can be approximated by Ae’-e*. But the steady 
state of x is e’, and Ax = } grad f(x), so that in the steady state 2A; = (grad f) -x. 
Since the components of grad f are available from the computation of x, \; can be 
obtained by performing the n multiplications necessary to form the dot product. 
The amount of equipment necessary for this is: 

n summing amplifiers, whose outputs are the components of grad f. 

n inverting amplifiers, whose outputs are the negative of the components 

of grad f. 

One summing amplifier to form the dot product. Multiplying equipment to 
compute the terms of the dot product x-(grad f). If multiplying servos 
were used above, one additional potentiometer per shaft is required. 

Two precautions must be taken in order to allow for the dynamic limitations 
of the multipliers. First, the problem may have to be slowed down in order that 
the multipliers remain within their rate limits. Second, the constant k may have 
to be adjusted carefully, in order to insure stability of the computing loop. 

Once an n X n matrix problem has been programmed, this programming—and 
the associated patching—can be used for any other n X n matrix, with only minor 
changes. Since the x; are restrained to lie between 1 and —1, the scaling need not 
be changed. Only the potentiometer settings and the sign of the variable being fed 
into the potentiometers must be adjusted. 

As an illustration, consider the following 3 X 3 matrix: 


. 2 2% 
Azi5 0 —3 x =| Ze 
2 aan) 1 3 


f(x) = Ax-x = 32° + 102422 + 4x23 — 6xe%3 + rs 
6%, + 10%. + 42; 


gradf=|107, — 62; 
47; _ 622 —_ 223 
Let k = 20. 


Equation (1) becomes 
% = 62, + 1022. + 423 + 2er; 
Z = 10x, — 623 + Were 
as = 4x7, — 622 + 223 + 20er;, 
where 
e = sgn (1 — (a + 22 + 25 )). 


The programming for this example is shown in Figure 1. The problem has been 
slowed down by a factor of twenty. 

The steady state of this system yields the eigenvector e' corresponding to the 
largest eigenvalue ), . 
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In order to compute the eigenvector e* corresponding to the smallest eigenvalue 
h; , equation 3 is used: 
(i, = —G6a, — 10x, — 423 + Wer 
14 = —102, + 6x; + 20er2 
(a; = —4a, + 6x2 — 2x3 + Wer. 


Finally, in order to compute the eigenvector e’ either of the two above systems 
of equations must be modified by the addition of an extra term. For example 
d = 62, + 10x. + 423 + 20er — fia 
Ze = 102; > 623 aa exe —_ foe, 
ds = 4x, — 622 + 223 + Wer; — a, 
where 
i 
f&}=e and « = sgn (2 & + 22k + 23 &). 
& 


5. Results. The eigenvalues and eigenvectors of two matrices were computed 
at Project Cyclone, using REACs and Reeves multiplying servos. One of the 
matrices was the 3 X 3 listed above. The other was a 6 X 6 matrix, listed in refer- 
ence 1 and reprinted below: 





.06667 .02634 — .04640 — .07368 — .02131 — .00431 | 
.02634 .26841 — .02243 .15952 — .05923 — .12797 
— .04640 — .02243 .10932 .05150 — .04100 08558 
— .07368 . 15952 .05150 25152 — .01141 — .07169 
— .02131 — .05923 — .04100 — .01141 14403 .01105 
| — .00431 — .12797 .08558 — .07169 01105 . 19450 | 





A comparison between the computed and actual values of the eigenvector 
components are listed in Table 1. 

The eigenvalues were not computed on the REAC, but a comparison may be 
made between (Ae‘-e’)/(e’-e') (where e’ is the computed eigenvector) and ),;. 
This is shown in Table 2. 


6. Additional Remarks. In the case of multiple eigenvalues, the eigenvectors 
are not, of course, unique. The method of this paper will then yield some set of 
orthonormal eigenvectors belonging to the eigenvalue. If two eigenvalues differ 
only slightly, the corresponding eigenvectors will be computed with less accuracy. 
But the corresponding eigenvalues will be obtained with no sacrifice in accuracy. 
In fact, the eigenvalues are relatively insensitive to errors in the eigenvectors. 
This is vividly demonstrated in Table 2, where the eigenvalues are seen to be cor- 
rect to five or six significant figures. 

The method of this paper can be extended to complex Hermitian matrices. 
The theorems listed in the mathematical formulation are basically still true if 
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TABLE 1 
Matrix 1 
e 











actual 


— .5194--- 
.7101--- 
AT5A--- 


actual 


.d261-- 

. 2639: - 

| .2134-- 
— .3044-- 
— .8006- - 
.2116-- 


actual 








| 


Actual 





computed actual computed actual computed 
. 7933 .7930--- | —.3189 | —.3184--- — .5193 
.6077 .6078--- | .3562 3555+ .7096 
— .0419 | —.0415--- | —.8794 — .8788--- .4762 
Matrix 2 
el e 
winnie’ a ms ual me canes actual computed 
.0408 .0414-- — .3558 — .3554: - .3253 
— .6762 | —.6763-- — .1227 — .1221-- . 2627 
.0387 .0394- - .6012 .6014- - .2129 
— .5746 | —.5745-- .5053 .5051- - — .3048 
. 1376 . 1382: - — ee | —.e-- | —.i4 
4361 -4361- - .4836 | -4838-- | .2104 
e e& 
gi Ts cial : leoeaiah E4 actual ‘enliaads 
— .38742 | —.3732-- .2977 . 2989 - . 7326 
— .6208 | —.5213 — .3411 — .3406- - — .2653 
.0686 | .0678- - — .6013 | —.6014-- F .4741 
.0214 .0222: - .5271 .5271-- . 2095 
—.4721 | —.4728-- —.2775 | —.2772-- .1771 
— .6004 | —.6002-- .2789 .2799-- — .3039 
TABLE 2 
Matrix 1 Matrix 2 
aes — a Se eee 
el ih. 7 —— 
——| iecielibmenall — 
Mi 6.727878 | 6.727881--- | M1 .4982334 
Ae 2.938089 | 2.938091--- Ae . 25946618 
Ms | —5.665964 | —5.665972--- | rs | .1759010 
| | Ne 0823481 
As | 01581329 
rs | 00268708 


| 
| 
| 
| 
| 
| 
| 


-4982344- - - 


.7328: > 
— .2650-- 
.4743-- 
- 2093 - - 
<eaaa 


— .3040- - 


. 25946632: - - 


. 1759013: - - 
-0823483 - - - 
| .01581310--- 


| .00268704- - - 





one extends the definition of inner product to complex vectors by the formula 


xy = Din wg [3). 


The quadratic form Ax-x = )>?j-; a;;0,4; is real, and the extremum property 
can be used in the same way to find the real and imaginary parts of the eigenvectors. 
Even normal matrices can be handled if one separates them into their real and 


imaginary parts [3]. 
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7. Conclusion. The eigenvalue problem for real symmetric, or Hermitian, 
matrices can be solved on an electronic analog computer by formulating it as an 
extremum problem. Both the eigenvectors and eigenvalues can be obtained. With 
care, three place accuracy can be obtained for the eigenvector. If the eigenvalue 
is computed by hand from the eigenvector, six place accuracy can be obtained for 
the eigenvalues. 


Reeves Instrument Corporation, 
Garden City, New York 
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Graphical Evaluation of a Convolution Integral 


By T. Mirsepassi 


1. Introduction. The analytical expression which defines the response of a 
linear system to an arbitrary excitation is, in general, in the form of a convolution 
integral [1, 2, 3], i.e., 


i) 2) = l Hie) -Alt = 9 o [ f(t — r)-h(r) dr = [ 10 - owe 


where: 

g(t) = response of system to unit step, and g(0) = 0, 

h(t) = g’(t) = response of system to unit impulse, 

f(t) = arbitrary excitation function, 

x(t) = response function of system to f(t), 

t = independent variable: time in the case of a time impulse such as in 
thermal [2] or electrical [3, 5] transients, and position in the case of space 
impulses, such as in deflection of beam or stretched cord under a space- 
variable load [4]. 

Usually the system—and therefore g(t) or h(t)—is known and it is desired to 
find the response to a given excitation. Sometimes the excitation function is not 
known and is to be found from a given response. An example of this type occurs 
in the conduction of heat when, in a given solid, the transient temperature at a 
point inside the solid is recorded and the temperature-time function on the boundary 
(or in the ambient) is of interest. Whenever analytical integration of (1) is possible, 
numerical evaluation can be carried out easily; otherwise it becomes unusually 
lengthy and uneconomical. For such cases graphical treatment is advisable. The 
classical method [3] of graphical integration of (1) consists of: 

1. Plotting h(7r), curve H in Fig. 1A 

2. Folding f(r), i.e., plotting it with +7 axis to the left, curve F in Fig. 1B 

3. Translating Fig. 1B on Fig. 1A such that the 7 axes coincide; for evaluating 
x(t) at t = t,, slide Fig. 1B on Fig. 1A until the origin of the 7 axis in Fig. 1B falls 
on rt = t, of the 7 axis in Fig. 1A, see Fig. 1C. 

4. Plotting the product curve, HF in Fig. 1C. Since HF represents f(t, — 7): 
h(7r), the value of x(t) for ¢ = t, is shown by the area under this curve and between 
the two lines 7 = 0 and r = ¢#,. The numerical value of this area, when carried 
along the ordinate of ¢ = 4, , locates one point of the response function, namely 
z(t,), P in Fig. 1D. Thus the mathematical process of convolution may be inter- 
preted graphically by four operations: folding Fig. 1B, translating Fig. 1C, then 
multiplying and integrating. The above method is tedious—for each additional 
evaluation of x(t) three operations (namely, translation, multiplication, and 
integration) must be repeated. 

In this paper a new method is described which is based on a finite-difference 
form of (1) and graphical multiplication. By means of this method, evaluation of 
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a(t) at any value of t is reduced to adding a series of readings. The method becomes 
especially time-saving when, on a given system, (1) is to be solved a number of 
times. 


2. Method. The function A(r) is plotted in Fig. 2A. In Fig. 2B, drawn 
on transparent paper, folding of the excitation f(r), drawn as a dotted line, is 
approximated by a step function drawn as a solid line; the duration of all steps along 
r axis is the same and is denoted by Ar, and the magnitude of the steps is denoted 
by F;, (¢ = 1, 2, --- ). The scales along the abscissa and ordinate in Fig. 2B have 
been taken, respectively, as equal to the scales along the abscissa and ordinate of 
Fig. 2A, and unit scale along the ordinate equals unit length. Now, superimpose 
Fig. 2B on Fig. 2A so that the horizontal axes coincide, see Fig. 2C, and that: 
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00’ = t = nAr, t, being a known value of ¢ at which the convolution integral is to 
be evaluated. Thus (1) may be written as follows: 


x( h = nAr) 


th nar 
| f(t — r)-h(r) dr = [ f(t — r)-h(r) dr 
0 0 


bo 
Il 


Ar 2Ar 
[ Fy-h(r) dr + | Fyaith(r) dr + --- 
0 Ar 


iar nr 
+ | F,-inath(r) dr + -+- + F,-h(r) dr. 
( 


i—1)Ar (n—1) Ar 
Since F; remains constant in each integral, it may be taken out of the integral sign; 
therefore 


7=n iAr =n 
(3) a(nAr) = > F.i41 / h(r) dr = = F,~i41 {g(iAr) — gl(t — 1)Ar]}. 
i=] (i—1) Ar i=l 


Let 
(4) G; = g(iAr) — gl(i — 1)Arz]. 


Then (3) becomes 


1=n 


(5) a(t, = nAr) = >> Py_ins-G;. 
i=l 
For graphical multiplication of F,,_;4:-G; , Fig. 2D is arranged. In this figure the 
center lines of Ar intervals are graduated by the scales 
(6) 8; = 8/G, 
where S = unit length. From these graduations one can read directly the product 


F,-:41°G; at the intersection P; , (Fig. 2D), of the step F,_;,; and the center line 
A;B;. In fact, since 





A;P; = F,-i41 x S 
therefore: 
(7) Reading at P; = A;P:/8; = (Fo—ia1-8)/(8/G,) = Friar Gi. 
Consequently (5) becomes 


=n i="n 


(8) a(t, = nAr) = >> Fyinn-G; = D> Reading at P;. 
t=] 


= 


Thus, by means of Fig. 2D, the evaluation of x(t) at 4; = nAr is reduced to the 
addition of n readings. 

The evaluation of x(t) for other values of ¢t is done similarly by sliding the trans- 
parent Fig. 2B on Fig. 2D to the new position and adding the readings at the new 
intersections of the steps and the center lines. 


3. Notes. 
1. The stepwise approximation of the excitation function may be omitted when- 
ever the function is approximately linear inside each interval. In fact, the average 
step, if drawn, will intersect the center line at about the same point as the function 
itself. The error due to stepwise approximation, however, will exist. 
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2. The response at t = nAr to unit step applied at ¢ = 0 is 


(9) g(nAr) = >> G; = >> Reading at B; 
t=] t=l 


where B; is the point of intersection at which the unit ordinate line intersects the 
center line of the 7**Ar interval, Fig. 2D. 


3. The response at t = (2 — })Ar to unit impulse applied at ¢ = 0, is 
(10) G;/Ar = (Reading at B;)/Ar. 


In fact, from (4) 
G;/Ar = {g(tAr) — g[(¢ — 1) Ar]}/Ar. 


The right side of this equation is the central difference of g’(¢) and thus represents: 


g(t) = h(t) 
at 


t = (4 — 4)Ar. 


Fig. 2D, which represents values proportional to the impulsive response, henceforth 
will be referred to as the “unit impulse chart”’. 

4. In problems concerning the conduction of heat, it is sometimes of interest 
to know the rate of temperature change at a point inside a solid when the surface 
(or ambient) temperature varies with time (e.g., quenching). This graphical 
technique can also be employed for problems of this type, i.e., problems where the 
derivatives of x(t) are of interest. In fact, using the Leibnitz rule [6]: 


z(t) = g f(r) -h(t — 7) dr |= “j(r) “W(t — rt) dr + f(t)-h(0) 
dt |_ Jo 0 


since g(t) represents the response at a known point inside the solid, which is initially 
at a steady state, when its surface (or ambient) temperature has undergone a unit 
step change—i.e., since: g’(t) = h (t) = O fort = 0; therefore 


(11) 2'(t) = l fits ~s} de = l f(t — r)-W(r) dr. 


In this case, (4) becomes 
(12) H; = h(tAr) — h{(t — 1)Az]; 


and, as explained before, a chart can be arranged similar to Fig. 2D. By a discussion 
analogous to that of Note 2, the resulting chart may be referred to as the “unit 
doublet chart”’. 

5. In complicated systems, g(t), h(t), or h’(t) may be obtained experimentally 
or, possibly, by an analog computer. Although, in such systems, response to an arbi- 
trary excitation may be found similarly, it often pays to make a unit-step run, 
arrange the “unit impulse chart”, and find x(t) graphically. This has two advan- 
tages: 

(a) Considerable saving of computer time 

(b) Possibility of solving additional response problems when the experimental 
setup or analog computer is no longer available. 
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4. Example. A linear system is defined by its response to the unit step according 


to the following relationship 
(13) g(t) =1—e". 
A. RESPONSE TO A GIVEN EXCITATION 
Find the response of this linear system to 
f(t) =, O0Ostsl 


f(t) = 1, 1 


IIA 


t<2 
1. Preparation of “unit impulse chart’’ 


Equation (4) becomes: 


G; = {1 - 


(14) e***) at {l —_— "sie _ g tie Be: gq. 


With due consideration for the required accuracy, select Ar for example, take 


Ar = 0.2 and calculate Table 1. 


TaBLeE 1 
Calculation of 1/G; from Equation (14). 





















































i 1 | 2 3 4 5 6 7 8 9 10 
iAr 0.2 | 0.4 | 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 
e~i-4r 1.000 0.8187; 0.670 | 0.549 | 0.4493 | 0.3679 | 0.3012 | 0.2466 | 0.2019 | 0.1653 
e~iAr 0.819 | 0.6703) 0.549 | 0.449 | 0.3679 | 0.3012 | 0.2466 | 0.2019 | 0.1653 | 0.1353 
G; 0.181 0.148 ; 0.121 | 0.100 | 0.0814 | 0.0667 | 0.0546 | 0.04470, 0.0366 | 0.0300 
1/G; 5.525 | 6.757 | 8.264 |10.000 {12.285 (14.992 (18.315 (22.371 (27.322 (33.333 
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TABLE 2 


Comparison of values of response obtained by evaluation of Equations (15) and (16), 
jinite-difference Equation (8), and actual readings of chart. 








j 


Ao oie | | 
‘ 0.2 0.4 0.6 0.8 j 1.0 1.2 1.4 | 1.6 1.8 2.0 
Response ‘NY | | 





Theoretical: Equations |0.0187\0.0703)0.1488 p.2409)0.2679 (0. 4824/0. 5763/0 .6531|0.7159|0.7675 
(15) & (16) | 
Finite-Diff. Equation (8) |0.0181/0.0692)0. 1473/0. 2475|0. 3658)0 . 4808|0. 5750/0 .6522|0.7155)0.7670 























Error (%) —3.2 |—1.5 |—1.0 |—0.7 |—0.6 |—0.3 |—0.2 |—0.1 |—0.05|—0.01 
Reading of Chart 0.018 |0.069 |0.147 |0.248 0.367 [0.480 10.575 (0.651 |0.714 |0.766 
Error (%) —3.7 |—1.8 |—1.2 |—0.5 |—0.2 |—0.5 |—0.2 |—0.3 |—0.2 |—0.2 
Graphical Error (%)* |—0.5 |—0.3 |—0.2 |+0.2 |+0.4 |—0.2 0 0.2 |—0.15)—0.19 
5 “Graphical error” = “Reading of Chart” error — “Finite Difference” error. ai 


Mark Ar intervals on the 7 axis in Fig. 3A; use an arbitrary scale. 
Select a unit length, S, (Fig. 3A), and graduate each center line of Ar intervals 


by the corresponding unit scale S;, where 
8; = (1/G,)8. 
For example, for 7 = 1 (the first center line), S, = (1/0.181)S = 5.5258, or 
0.18, = 0.5525-8. 
2. Evaluation of response 


Plot the folding of f(t), take the same Ar intervals as in Fig. 3A, and use unit 
scale § for the ordinate; see Fig. 3B which should be drawn on transparent paper. 

Superimpose Fig. 3B on Fig. 3A, as is shown by the dotted line in Fig. 3A for 
t = 0.8, and add the readings; thus, for ¢ = 0.8, the response is: 


x(t) = 0.010 + 0.036 + 0.076 + 0.127 = 0.248. 


Repeat this second step for other values of ¢ and plot Fig. 3C. It is interesting to 
analyze the errors in this example. The analytical solution is: 


t 
(15) tt) = [ soy aie ¢4 5" 44, be OO 


lA 
IIA 


t-—1 t 
(16) 2(t) = | l-e dr + / (t—r)e' dr =1—(e—1)e’, for t= 1. 
0 t-1 


Table 2 shows theoretical values computed from (15) and (16), finite difference 
values computed from (8), and values obtained by reading the chart. 
It is interesting to note that the error is primarily introduced by the difference 
approximation, and that the graphical error does not exceed 0.5% . 


B. EXCITATION FUNCTION FOR A GIVEN RESPONSE 


Find the excitation function, f(t), if in the above system a response, x(t), is 
given as follows: 








| t | 0 0.2 0.4 0.6 0.8 1.0 1.2 4 1.6 1.8 22-0 
| 


x(t) | 0 10.0180. 069/0. 147/0.24810.36710.48010.575)0.651/0.714|0.766 
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Evaluation of excitation function 


Plot x(7r) with the 7 axis increasing from right to left with the abscissa scale the 
same as in Fig. 3A, (i.e., the distance between two center lines representing Ar = 
0.2), and the ordinate scale equal to the unit length § (Fig. 4A, which should be 
drawn on transparent paper). 


Superimpose Fig. 4A on Fig. 3A as shown in enlarged view in Fig. 4B. Let A 
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Figure 4D 


denote the intersection of the unknown f(t) with the center line. From (5) written 
forn = 1, 
x(Ar) = F,-G, . 


The lefthand side of this equation is the reading at P from the graduated unit 
scale. The righthand side is the reading at the point A from the graduated center 
line. The lefthand side is known—e.g., 0.018 in Fig. 4B. The righthand side must 
then equal the reading at P. Therefore, the location of A is established in this case 
as the point 0.018 (read from the graduations of the center line) along the center 
line. 

Then, slide Fig. 4A on Fig. 3A to the second position, as shown in Fig. 4C. Let 
B denote the new intersection of f(t) and the first center line. From (5), written 
forn = 2: 

x(t) = F.-G, + F,-Ge ° 


Since: x(t) = Reading at Q from graduated unit scale (i.e., 0.069 in Fig. 4C) and 
F,-G, = Reading at A from graduated center line (i.e., 0.0148 in Fig. 4C), 
therefore: 


F,-G, = Reading at B from the graduated center line 


x(t) — F,-G, = 0.0690 — 0.0148 = 0.0542 


and in Fig. 4C, Point B is located where the reading along the center line is 0.0542. 
Thus, after each translation, the reading at the new intersection of F(t) and the 
first center line is obtained by subtracting from z(t) (read from graduated unit 
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scale) the sum of the readings at A, B, --- etc. The complete result is shown in 
Fig. 4D. 


5. Concluding Remarks. The method described in this paper is based on a 
finite-difference form of the convolution integral and a graphical multiplication. 
Once a chart is arranged, an evaluation of a convolution integral reduces to adding 
a series of readings obtained from the intersections of a line and a number of gradu- 
ated scales. The method is especially time-saving when, on a given system, the 
convolution integral is to be solved many times. Heat-transfer charts based on this 
method are under preparation; the first part appeared in [7] and the continuation 
will be published as completed. 
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University of California Extension, 
Los Angeles, California 


1. W. T. Tomson, Laplace Transformation, Prentice Hall, Inc., New York, 1950, p. 37-38. 

2. H.S. Carstaw & J. C. JazGer, Conduction of Heat in Solids, Oxford Press, New York, 
1950, p. 20. 

3. M. F. Garpner & J. L. Barnes, Transients in Linear Systems, John Wiley & Sons, 
Inc., New York, 1942, p. 231-234. 

4. L. A. Pipes, Applied Mathematics for Engineers and Physicists, McGraw-Hill Book Co., 
New York, 1946, p. 213 and 222. 

5. Stanrorp GoLpMAN, T'ransformation Calculus and Electrical Transients, Prentice Hall, 
Inc., New York, 1950. 

6. C. R. Write, Jr., Advanced Engineering Mathematics, McGraw-Hill Book Co., Inc., 
New York, 1951 p. 591. 

7. T. J. Mrrsepasst, ‘‘Heat-transfer charts for time-variable boundary conditions—semi- 
infinite solid,’’ British Chemical Engineering, March 1959, p. 130-136. 











wn in 


on a 
ation. 
dding 
radu- 
1, the 
n this 
ation 


37-38. 
York, 


Sons, 
1k Co., 
e Hall, 
, Inc., 


—semi- 





The Use of the Central Limit Theorem for 
Interpolating in Tables of Probability 
Distribution Functions 


By Gerard Salton 


In using tables of probability density and distribution functions, the difficulty 
of interpolating for functional values which are not directly tabulated constitutes a 
major problem. In particular, when the variance of a random variable becomes 
small, the corresponding density and distribution functions approach, respectively, 
a delta-function and a step-function. As a result, for small variations in the vari- 
ables, there occur very large variations in the corresponding functional values, and 
interpolation by difference methods will give very poor approximations. Moreover, 
space limitations in a volume of tables often make jit impractical to tabulate a given 
function on a mesh fine enough to allow for accurate interpolation. This often 
increases the difficulty of interpolating for a given probability distribution. 

By a fundamental limit theorem of probability theory it is known that, under 
rather general conditions, the sum of a set of n independent random variables 
appropriately standardized by the mean and the standard deviation is asymptoti- 
cally normal with mean 0 and variance 1, as n tends to infinity [1]. A method which 
makes use of this property to improve the accuracy of interpolating in tables of 
probability distribution functions was suggested by Rossow [2]. The binomial 
probability distribution is chosen here for purposes of illustration. The method 
described is however applicable to any set of random variables which obeys the 
central limit theorem. 

The cumulative binomial probability distribution may be denoted by 


E(n,r, p) = 2 C,"(1 — p)” “p' 


where 0 S p S landO <r & n. If one considers a series of n independent repeti- 
tions of some random experiment, then E(n, r, p) represents the probability of at 
least r successes in n repetitions of the experiment, p being the probability of success 
for each experiment. By the theorem previously stated, it is known that for a given 
value of p, and for increasing values of n 


1 
V2" 
where x = (r — np — 3)/+/np(1 — p), and N(x) denotes the cumulative standard 
normal distribution. The effectiveness of the approximation increases with increas- 
ing values of +/np(1 — p), so that for a given value of n, the approximation is 
closest for p = 34. The argument zx, corresponding to a given value of N(z), is 
often called the normal deviate corresponding to the total frequency N. 

Consider now the problem of interpolating in a table of the binomial probability 


E(n,r, p) © / edt = 1 — N(z) 
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distribution. As an example, let it be desired to find some value E(n , ro , Po + Ap), 
where the function is not tabulated for the argument p) + Ap. The standard 
procedure consists in taking the values of E(n, r, p) corresponding to the arguments 
NoroPo , MooP1 , *** » MoPe+1 and in constructing a k + 1** degree polynomial which 
takes on the given functional values at the k + 2 given points. The value of the 
polynomial at the point noropo + Ap is then calculated and taken to be equal to 
E(m , To, Po + Ap). 

It is proposed to modify this procedure by making use of the fact that the 
binomial distribution is often much better approximated by the corresponding 
normal distribution than by the interpolation polynomial which is constructed in 
the standard method. To this effect, it is first assumed that E(m , 70, po) = 1 — 
N(x), E(mo, m%, Pi) = 1 — N(x), --+, E(mo, m0, Perr) = 1 — N(ae41). The 
k + 1 normal deviates corresponding to the k + 1 values of 1 — N are then found 
in a table of the normal distribution. Thereafter the interpolation is performed in 
the normal deviates, that is, from the values —2x; corresponding to the arguments 
pi, a value —2z, is determined which corresponds to the argument pp + Ap. The 
value of 1 — N(z,) is then looked up in a table of the normal distribution, and 
E(m , tT, Po + Ap) is approximated by 1 — N(x,).* The procedure is outlined in 
the following diagram, where the arrows denote a table look-up and the brace 
stands for the evaluation of an ordinary interpolation polynomial: 


po > E(m,%m,~) = 1— N(m) — —%» | 
1 — N(x) aS 7a 


pr — E(nm, 1,7) 
—2x(po + Ap) 1 — N(z,) 


Piri — E(m, 10, Presi) = 1 — N(aey1) > —Teqr - 
= E(m, 10, po + Ap). 


The procedure is easier to perform if tables are available which tabulate the 
normal deviates as a function of the cumulative normal probabilities [3]. However, 
tables which give N(x) as a function of x can also be used [4]. In either case, the 
accuracy of the method depends on the availability of accurate values of the normal 
deviates; this is especially important for cumulative probabilities close to 0 and 
close to 1 where the deviates change rapidly for small changes in the probabilities. 
It should be noted that the additional labor required by the proposed method is 
restricted to some table look-up operations which can be performed rapidly. 

A study was made of the efficiency of the method, based on a recent tabulation 
of the cumulative binomial probability distribution [5]. It is stated in that volume 
that interpolation by standard difference methods gives poor results for certain 
ranges of the arguments. In particular, p-wise interpolation is poor forn > 100 and 
p small, r-wise interpolation is poor for n < 100 and p small, and n-wise interpola- 
tion is inaccurate for n > 100 and p close to 3. 

The normal deviate method fills the need for both p-wise and n-wise interpola- 





* Since the arguments of 1 — N(z) are the negatives of the corresponding arguments of 
N (a), the method can also be carried out by interpolating in the values x; , rather than in the 
values —z; . Such an interpolation will result in the determination of the function N(z,) cor- 
responding to 1 — E(mo, ro, po + Ap). 
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tion, since it is most accurate for large n where the standard method fails. It can 
also be used for interpolation in r when n is not too small. When interpolating for 
p and for n, a four-point interpolation in the normal deviates was found to result in 
final probabilities whose error did not exceed 5.10°°, except for very small p 
(p < 0.03) where the error could be as high as 5.10~*. Four-point interpolation in 
r for n > 30 resulted in errors not exceeding 5.10“. 

A different method which also improves the interpolation in tables of prob- 
ability distribution functions consists in approximating the given distribution by 
the logistic function [6] 

1 


Mz) = Te 


in lieu of the cumulative standard normal distribution. In fact it is known that 
L(x)  N(z2), 


the two functions differing by less than 0.01 when 8 is close to 1.7. The interpolation 
proceeds exactly as before except that L(x) is used for N(x). The arguments z 
corresponding to various values of L(x) are sometimes called logits. By inverting 
the expression for L(x) it is seen that 


t= di a#.. . 
B 1-L 
Tables of the natural logarithms and exponential tables can therefore be used, if 
logit tables are not available (6). 

Linear interpolation in the logits will sometimes give slightly better values than 
the corresponding interpolation in the normal deviates. However, a higher order 
interpolation formula brings much less improvement in the logits than it does in 
the normal deviates. A four-point interpolation in the normal deviates is usually 
to be preferred to the corresponding four-point formula in logits. It is again im- 
portant to use accurate logit tables, especially for values of L close to 0 and 1. Since 
tables of the logistic function may not be generally available, interpolation in the 
normal deviates is usually easier to perform and will moreover give more accurate 
results for higher order interpolation formulas. 

The following examples give a good idea of the improvement which can be 


Comparison of Interpolation Methods 




















E (600,225,0.375) E (1000,53,0.0625) E (39,6,0.08) r-wise, | E (689,285,0.42) 
p-wise, Ap=0.5 p-wise, Ap=0.25 Ar=1 | #-wise, An=39 
Correct Values .51541 .90685 -08786 .64612 
Ordinary 2-pt. .51475 (—66) | .87776 (—2909) | .11597 (+2811) | .63732 (—880) 
Ordinary 4-pt. .51531 (—10) | .90178 (—507) -09194 (+408) .64157 (—455) 
Normal ” ' 
Deviates 2-pt. .51524 (—17) | .90497 (—188) .09004 (+218) -64263 (—349) 
Normal ” 
Deviates 4-pt. .51542 (+1) -90685 (0) .08770 (—16) .64616 (+4) 
Logits 2-pt. .51537 (—4) .91249 (+564) .08456 (—330) -64401 (—211) 


Logits 4-pt. -51543 (+2) -90973 (+288) | .08726 (—60) -64870 (+258) 
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achieved over the standard interpolation method by using either normal deviates or 
logits. In each case, the error in the value obtained is shown in parentheses. 


Computation Laboratory, Harvard University, 
Cambridge, Massachusetts 


i. Haratp Cramér, Mathematical Methods of Statistics, Princeton University Press, 
Princeton, 1951, p. 213-220. 

2. E. Rossow, Technische Universitit, Berlin, private communication. 

3. R. A. Fisoer & F. Yates, Statistical Tables for Biological, Agricultural and Medical 
Research, Oliver and Boyd Ltd., Edinburgh, 1948. 

4. Harvarp University, CompuTaTion LABoraTory, Annals, v.23: Tables of the Error 
Function and of its First Twenty Derivatives, Harvard University Press, Cambridge, Mass., 
1952. 

5. Harvarp UNIvEersity, CompuTaTION LasporaTory, Annals, v. 35: Tables of the Cumu- 
lative Binomial Distribution, Harvard University Press, Cambridge, Mass., 1955. 

6. JosePpH Berkson, ‘‘A statistically precise and relatively simple method of estimating 
the bio-assay with quantal response, based on the logistic function,’’ Journal of the American 
Statistical Association, v. 48, No. 263, September, 1953. 
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34[A, B, C, D].—L. Fuavien, Nowvelles Tables Numériques pour les Fonctions 
Usuelles de ’ Analyse, Gauthier-Villars, Paris, 1958, 63 p., 21 em. Price $0.85. 
This little volume has been compiled for students who are preparing “aux 

grands écoles scientifiques.” It has been edited to conform with the 1956 program 

published by the ministére de l’Education Nationale. Tables included are: 
I. Avertissement. 

II. Values of n’, n*, and n™“, n™*, n™*, +/n V/n to 4D, for n = 1(1)1000. 

III. Mantissa of logy z to 4D for x = 1(1)999. 

IV. Log, x to 4D for x = 1(1)1000. 

V. 10° to 4D for x = 0(.01)1. 

VI. Sin z, tan 2, ctn x, cos x to 4D, in degrees for x = 0(0°.1)90°. 

VIL. Sin z, tan 2, etn x, cos x to 4D, in grades for x = 0(0°.1) 100’. 

VIII. Are sin x in radians to 5D for x = 0(0.01)1. Are tan z in radians to 5D for 

x = 0(0.01)1 (var.)100. 

IX. Radians r: to degrees (decimal), to grades (decimal), to 4D for r = 1(1)10, 
Radians to degrees, minutes, and seconds for r = 0.1(0.1)0.9, 0.01(0.01)0.09. 
0.001 (0.001 )0.009, 0.0001 (0.0001 )0.0009. 

X. Degrees d, minutes m, seconds s: to grades (4D), to radians (5D) for d = 
1(1°)90°, m = 1(1)60’, s = 1(1”)60”. 

XI. Grades G: to radians (4D), to (decimal) degrees (1D), to degrees, minutes, 
for G = 1(1°)100°. Centigrades C to minutes, seconds, (1D), for C = 
1(1°)100°. 

XII. “Remarkable” numbers (e.g., 7, e, n!) and their common logarithms. 


The Tables are well designed and easy to read. 
Ricuarp 8. BurIncTon 
Bureau of Ordnance, 
Navy Department, 
Washington, District of Columbia 


35[A, B, I, N]|.—P. Montaene, Tables abrégées de puissances entiéres, Dunod, Paris, 

1958, xv + 411 p. + loose appendix 32 p., 28 em. Price 5600 francs. 

These extensive tables are designed to facilitate the evaluation of x", where n 
is a positive integer, more accurately than is possible with eight-place logarithmic 
tables. The main tables are divided into three sections comprising tables referred 
to as small, medium and large. The values of x” are given in the small and medium 
tables to 15S, and in the large tables to 10-118. Numerous special signs attached 
to the last figure give some information about the next place. There are no differ- 
ences. The values of the arguments are: 

Small tables (a), pages 3-17: x = .2(.1)2, n = 0(1)m, where m depends on 

x in the following way: 


x 2 3 4 .5(.1).9 1.1(.1)2 
m 600 400 250 200 150 


Small tables (b), pages 18-32: « = .1(.01)1.53, n = 2(1)78. 
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Medium tables, pages 64-157: x = .097(.001)1.263, n = 2(1)26. 

Large tables, pages 160-411: x = 0(.0001)1.2603, n = 2(1)10. 
In all tables except small (a), x is the vertical argument, and there is an overlap 
of about 7 lines between pages, instead of the usual one line or even none. This 
accounts for initial and terminal arguments ending in 7 and 3; thus, the medium 
tables are basically for x = .1(.001)1.26, the three extra arguments at each end 
being added for convenience in interpolation and so forth. The tables are photo- 
graphically reproduced from manuscript written by professional calligraphers. 

A loose appendix contains 29 pages of Tables annexes and 3 pages reserved for 
notes. The varied contents defy brief description in full. About half the pages con- 
tain matter relating primarily to the main tables, such as illustrative diagrams, in- 
dexes, and tables indicating the number of differences needed in interpolation. The 
remaining pages contain matter of more general applicability. Outstanding are two 
tables of Everett interpolation coefficients, one giving exact coefficients of second 
and fourth differences at interval 0.002, and the other, phenomenally extensive, 
giving at interval 0.01 coefficients to 12-13D of all even differences up to the six- 
teenth; except in the case of coefficients of the 14th and 16th differences, even 
differences (up to the second or fourth) of the interpolation coefficients are also 
given. There are also rather extensive tables of factorials, reciprocals of factorials, 
binomial coefficients, and useful constants (including Bernoulli and Euler numbers). 

There are few references; it seems a pity not to mention the important British 
Association tables of powers [1]. The tables of x” will certainly be found useful for 
tabular arguments, but when much interpolation would be needed, the use of log- 
arithms as an alternative may still appeal on occasion. Obviously much work has 
gone into checking the tables, but a copious page of errata is pasted in, and it is 
stated that the last check revealed less than one error pér ten pages. There do exist, 
however, radix and similar special logarithmic tables which are slim, simple, and of 
absolute accuracy. Nevertheless, anyone much concerned with integral powers 
should consider what possibilities this volume may have for his purpose. An excellent 
feature of that part of the large tables which relates to the range 0 S x & 1 is that, 
for given 2, all the powers of x and 1 — z are contained on one line running right 
across two pages visible at an opening; this facilitates the solution of equations of 
the form 2*(1 — z)*® = k, which occur in connection with chemical equilibria. 
The extensions for z > 1 may be found useful in connection with problems on com- 
pound interest. 

Editions in English and several other languages are in preparation, as also is a 
table of fourth powers. 

ALF. 


1. British Association for the Advancement of Science, Mathematical Tables, vol. IX: 
Table of Powers giving Integral Powers of Integers, Cambridge University Press, 1940. See mrac, 
Review 169, v. 1, 1945, p. 355-356. 

36[G].—L. Lunetur & M. Scz, ‘Sulla ricerca dei k-archi completi mediante una 
calcolatrice elettronics,” Atti Convegno Intern. Reticoli e Geometrie Proiettive, 
Palermo-Messina, 1957 (Roma 1958), p. 81-86. 

-archi completi nei piant proiettivi desarguesiani 

di rango 8 e 16, Politecnico di Milano, Centro di Calcolo Numerici, Milano 1958, 

11 p. 
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The authors have programmed a CRC 102 A/P to search for complete k-arcs 
in a finite projective plane coordinatized by a Galois field. A k-are is a set of k points 
such that no three lie on a line. It is complete if it is not contained in a (k + 1)-are. 
To eliminate the most obvious duplications, only k-ares through the four funda- 
mental points (0, 0, 1), (0, 1,0), (1, 0,0), (1, 1, 1) were considered, but the authors 
regard as distinct k-ares those which result from one another under collineations 
permuting the four points. 

The procedures in the first paper are limited to planes of prime order. A descrip- 
tion of the program is given and the results of a complete search of the plane of 
order seven are announced. Examples of complete k-ares are given for planes of 
orders 11, 13, and 17. 

In the case of the plane of order seven, 40 6-ares are tabulated. 

In the second paper, the program is extended to utilize the irreducible poly- 
nomials z° + x + 1 and x* + x + 1 to calculate in terms of the coordinatization by 
GF(8) and GF(16) for planes of these orders. (A misapprehension as to the general 
suitability of polynomials x” + zx + 1 does not affect the results in these cases.) 
Ten 10-ares are listed for the plane of order eight, and 45 6-ares. If permutations 
were considered, the authors could have reduced the tabulation to three 10-arcs 
and fifteen 6-arcs. Indeed, as they point out, the case of the 6-ares can be reduced 
to the statement: any 5-are can be completed to a complete 6-are in precisely three 
ways. For the case 16 no exhaustive tabulation is given. Some 18 10-ares, one 11- 
are, two 12-ares and three 18-ares are listed. The last do not contain conies (ovals) ; 
this answers a question of B. Segre as to the existence of such ares. 

J. D. Swirr 


University of California, 
Los Angeles, California 


37[H, I, S]|.—R. L. Coampers & E. V. Somers, Solution of a radiation type non- 
linear differential equation, Scientific Paper 8-0529-P6, Westinghouse Research 
Laboratories, Pittsburgh, Pennsylvania, October 1958. 


The treatment of a heat transfer problem [1, 2] dealing with radiation from a 
cooling fin gives rise to a second order non-linear differential equation 


d’°0 1 dé 

tsp pe * 
with boundary conditions @ = 0 at r = 0, d@/dr = 0 at r = 1. Tabulated results 
are presented from the numerical solution of the above equation over a range of 
values of r from 0 to 1.0, for the parametric values of p = 1.001, 1.1, 1.25, 1.5, 2.0, 
and 3.0 and y = 0 to 4. The tables were computed on an IBM 704 digital computer 
by a finite-difference method. A copy has been deposited in the UMT file. 


= 0; 


H.P. 


1. Rosert L. CuamBers, The determination of radiation fin efficiency and temperature dis- 
tribution for one-dimensional heat flow in a circular fin, University of Pittsburgh thesis, 1958. 

2. Ropert L. CuampBers & E. V. Somers, Radiation fin efficiency for one-dimensional heat 
flow in a circular fin, Westinghouse Scientific Paper 8-0529-P4, 1958. 


38{I, X, Z|.—R. M. Pearce, Digital Computer Solution of the Two Group Diffusion 
Equations in Cartesian or Cylindrical Geometry with Application to the Datatron, 
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Report AECL No. 487, 1957, 27 p., 27 ¢.m. Available from Scientific Document 
Distribution Office, Atomic Energy of Canada Limited, Chalk River, Ontario, 
Canada. Price $1.00. 


This report derives a set of difference equations for use in solving the two-group 
diffusion equations in XY or RZ geometries, and describes a Datatron program to 
solve the difference equations. The Datatron program and the output from a sam- 
ple problem are included. 

The Datatron program can handle problems with a maximum of 324 mesh points. 
This is less than some other currently available programs can handle. 

The boundaries may be black (the value of the flux is forced to zero) lines of 
symmetry, or what is referred to in the report as reflector termination. In the 
case of reflector termination, an attempt is made to simulate the effects of a par- 
tially reflecting material bounding the cell. 

The iteration scheme uses point relaxation. The latest values are used except in 
one case where their use would require extra machine time per iteration. No attempt 
is made to accelerate the convergence of the iteration scheme, even though there 
are several methods available which have been used successfully in similar pro- 
grams. 

CHARLES DAWSON 
Applied Mathematics Laboratory, 
David Taylor Model Basin, 
Washington 7, District of Columbia 


39[L].—H. V. McInrosu, A. Kueppner & D. F. Mrnnur, Tables of the Herglotz 
Polynomials of Orders 3-$ Transformation Coefficients for Spherical Harmonics, 
Ballistic Research Laboratories, Memorandum Report No. 1097, Aberdeen 
Proving Ground, Maryland, 1957, 162 p., 28 cm. 


The Herglotz polynomials Hj, (R) are the matrix elements of the unitary 
irreducible representations of the three-dimensional rotation group and are defined 
by: 


H".(o0, 6, y) = (—1)*“*e"e™*(cos ¢/2)" *™ (sin ¢/2)** 





. a (—1)* tan” $/2 = 
7 (n+k—s)'"(7 —k+8)'(n —j — 8)! 8!’ 
where a, ¢, y are the Euler angles of the rotation R. These polynomials have the 
useful property that if R takes the coordinate system r, 6, — into r, 6’, &’, then 


P,3 (cos 6’)e* = >> Hj.(R)P,* (cos 6)e™* 
k=—n 
where P,,’ (cos 0) e’* are the normalized spherical harmonics. 
In this table the values of Hj, (0, ¢, 0) are given for n = $(4)$ and@ = 0° 
(1°) 90°. 
AvuTHORS’ SUMMARY 
40(I, M, P|.—-R. L. Murray & L. A. Minx, Tables of " for Reactor Slabs, Cylinders, 
and Spheres, n = 1 ton = 20, Department of Engineering Research, North 
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Carolina State College, Raleigh, North Carolina, Bulletin No. 70, 1958, 56 p., 
28 em. Price $1.50. 


This paper contains a tabulation of average neutron fluxes for slab, cylindrical 
and spherical reactors. Specificaily, the following functions are tabulated to 8 
decimal places: 


éx/2 
Slab o"(z) = TB l [Cos z]" dx with zx = xz/H 
/ 


I 


Ww 550 
Cylinder ¢"(r) . of [Jo(x)|"2? dx with x = jor/Rs 
(8}0)? Jo 








‘ bx Fa: n 
Sphere ¢"(r) = ia l ‘= *| a dx with x = r/R. 
R and H represent the extrapolated boundries of the reactor. 

5, the fraction of the total dimension over which the average is performed, is 
given from 0 to 1 in steps of .01, and m assumes integer values from 1 to 20. The 
above values of ¢" actually represent o"/¢.", where ¢.” is the central flux normalized 
to unity. 

These tables were calculated on an IBM 650. Details of the methods used are 
presented, as well as a method for interpolating for values not tabulated. 

The paper illustrates the application of these flux averages by means of a few 
worked examples. 

Rosert Bropsky 
Department of the Navy, 
Washington 25, District of Columbia 


41[P, T, Z]|.—Ernest F. Jounson, “Automatic Process Control’—-Chapter II, 
Advances in Chemical Engineering, Thomas B. Drew & John W. Hoopes, Jr., 
Editors, Academic Press Inc., New York, 1958, x + 338 p., 23 cm. Price $9.50. 


This chapter should be a convenient primer for those to whom this subject is 
new. Terms are defined clearly, equations given and the bibliography is excellent. 
Some of the more advanced notions such as “three-mode control’ are discussed 
adequately. It is unfortunate that the topic of “sampled-data systems” is no more 
than mentioned, as this is an area of equal importance functionally, and will be of 
much greater importance in the future when chemical engineers finally appreciate 
that digital computers are much more powerful and flexible than analog systems. 


GILBERT W. KING 
IBM Research Center 
Yorktown Heights, New York 


42(W, X, Z]|.—Armour Researcu Founpation, Proceedings of the Fifth Annual 
Computer Applications Symposium, 1958, Sponsored by the Armour Research 
Foundation of Illinois Institute of Technology, x + 153 p., 23 em. Price $3.00. 


This small volume contains copies of the papers presented at the 1958 Computer 
Applications Symposium held on 29 October in Chicago and sponsored by the 
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Armour Research Foundation. A list of papers presented ai this meeting is given 
below. 

1. Operations Research and the Automation of Banking Procedures—R. A. 
BYERLY 

2. Information Systems Modernization in the Air Materiel Command—D. E. 
ELLETT 

3. Utilization of Computers for Information Retrieval—A. OpLER 

4. Problems and Prospects of Data-processing for Defense—C. A. PHILLIPS 

5. An Integrated Data-processing System with Remote Input and Output— 
R. D. WHIsLER 

6. The Role of Character-Recognition Devices in Data-processing Systems- 
R. L. HARRELL 

7. Input-Output—Key or Bottleneck?—R. D. ELBourN 

8. Scientific Uses of a Medium-Scale Computer with Extensive Accessory Fe: 
tures—R. A. HAERTLE 

9. The Design of Optimum Systems—R. R. Brown 

10. Computer Applications in the Numerical Control of Machine Tools 
R. B. Cieae 

11. Frontiers in Computer Technology—R. W. Hammina 

12. Computer Sharing by a Group of Consulting Engineering Firms—E. M. 
CHASTAIN AND J. C. McCatui 

13. Current Developments in Computer Programming Techniques—I’. Way, 
Ill 

14. The Future of Automatic Programming—W. F. Bauer 


ae A 


43[X].—Vera Ritey & Saut I. Gass, Linear Programming and Associated Tech- 
niques, Bibliographic Reference Series No. 5, Published for the Operations Re- 
search Office, The Johns Hopkins University, The Johns Hopkins Press, Balti- 
more, 1958, x + 613 p., 23 em. Price $6.00. 


This bibliography, a revised edition of BRS-2, Programming for Policy Decision, 
March 1954, includes over 1000 abstracts of articles, books, monographs, theses, 
conference proceedings, etc., dealing with linear programming and related topics 
such as nonlinear programming, dynamic programming, and game theory. The 
comprehensive work embraces references on the history, progress, and application 
of mathematical programming. It is divided into four sections: Part I, ‘‘Introduc- 
tion,” covers the early development and basic concepts of linear programming; 
Part II, “General Theory,’ embraces a wide range of topics such as advanced 
mathematical aspects of linear programming, computational methods and machine 
techniques, linear inequalities and convex sets, and theory of games; Part III, 
“Applications,” presents applications to industrial, agricultural, and military prob- 
lems, and contains a basic bibliography of material related to the field of production 
scheduling and inventory control; Part IV, “Nonlinear and Dynamic Program- 
ming,”’ covers the mathematical and computational aspects of nonlinear and dy- 
namic programming. Each section is prefaced by an expository discussion of the 
scope and contents of the material listed. 

Notwithstanding the spate of impressive publications on the applications of 





line: 
adv 
mal 
met 
solv 
tria 
pro; 
req) 


and 
and 
of | 


van 


Apr 
Dax 
Was 


pre 
lov 


wh 


ea- 


och- 
Re- 
Iti- 


ion, 
ses 
ies 
Che 
jon 
uc- 
ng; 
ced 
\ine 
[I], 
‘ob- 
lon 
4m- 
dy- 
the 


: of 





REVIEWS AND DESCRIPTIONS OF TABLES AND BOOKS 223 


linear programming, its record of accomplishment does not support the glowing 
advocacy of it as a practical tool in management science; its contribution has been 
marginal at best. The lack of efficient computational techniques and machine 
methods precludes the use of mathematical programming as an effective vehicle in 
solving the multifarious and vexing problems inherent in many large-scale indus- 
trial, governmental, and military operations. The reviewer does not discount the 
progress already made, but wishes to stress the fact that extensive research is still 
required to effectively solve the complex and formidable problems of management. 

In the opinion of the reviewer, the bibliography will serve as an invaluable aid 
and guide for obtaining a detailed account of the current mathematical techniques 
and applications. The broad spectrum of references describing the work in the field 
of linear programming makes the subject accessible even to readers without ad- 
vanced training in mathematics. 

Mitton SIEGEL 

Applied Mathematics Laboratory, 
David Taylor Model Basin, 
Washington 7, District of Columbia 


44(X, Z|.—_R. W. Merzcer, Elementary Mathematical Programming, John Wiley & 
Sons, Inc., New York, 1958, ix + 246 p., 22.5 em. Price $5.95. 


This book was written as an attempt to fill a gap which exists in the literature 
concerning mathematical programming. The author tries with singular success to 
hit the middle ground between purely literary discussions for the layman, which 
are of no technical value, and highly technical papers, which are incomprehensible 
to all except professional mathematicians. Mr. Metzger assumes his reader has a 
limited background in mathematics but wishes to understand the basic techniques 
and applications of mathematical programming. 

In Chapters 2, 3, and 4 Mr. Metzger gives an excellent exposition of the dis- 
tribution, simplex and approximation methods. The techniques are worked out by 
use of sample problems which are well chosen for their simplicity and applicability. 
The mathematics of the analysis and solution is complete, accurate, and under- 
standable. Mr. Metzger goes to considerable lengths to keep the reader from getting 
lost in unfamiliar mathematics. At no time is the reader informed that “it is not 
difficult to show that . . .”” Each step is carefully explained. No proofs are given for 
any of the methods. However, intuitive explanations keep the reader from feeling 
that the mathematics involved is a kind of “black box” which magically produces 
the right answer. Adequate references are given for those who wish to verify the 
proofs. Each method of solution is summarized in a step by step listing of the pro- 
cedure for easy reference. 

Chapter 6 is a short discussion of computers and their applications in this field. 
Existing programs for various computers are mentioned and some indications given 
as to size limitations and time required for computation. 

The remaining chapters are concerned with various applications of mathematical 
programming to industrial and business problems. The analysis preceding and fol- 
lowing the mathematical computation are emphasized and illustrated with some- 
what simplified but practical problems. 

The book is well documented, and contains a fairly good bibliography. It is best 
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suited for individual study but could be used as a text or supplementary text in a 
course. The only fault lies in the scarcity of practice problems other than those 
worked out in full in the text. 

This book will be extremely valuable to any person who wishes an introduction 
to mathematical programming, whether or not he is a mathematician. Mr. Metzger 
should be commended for fulfilling a real need. 


JEAN PoRTER 
9308 48th Avenue, 
College Park, Maryland 


45[Z]._-_Nep Cuapin, An Introduction to Automatic Computers, D. Van Nostrand 
Company, Inc., Princeton, New Jersey, 1957, viii + 525 p., 25 em. Price $8.75. 


The stated objective of this revised edition of a 1955 book is to present com- 
puters from the business systems point of view. This is developed through two sub- 
books, one set of chapters dealing with the logical and engineering nature of a 
computer and with structural and operating characteristics of a large number of 
commercial devices and computers, and another set of chapters dealing with the 
applications of computers in business operations. The programming chapters are 
related to computers, not to application. 

The style and depth are largely that of the numerous popularizations in business, 
systems and accounting periodicals. There is prevalence of many of the common 
clichés and generalizations, preoccupation with explaining or destroying ‘‘popular 
misconceptions”, and some narrowness of viewpoint about the nature of business 
management and the role of computers in it. Many managers are still seeking con- 
ceptual simplification and education about computers and their relation to manage- 
ment. It is doubtful that they approach this search as less than serious students, 
students hoping to find scientific classification and orderly development rather than 
demeaning and patronizing popularization. 

To its credit, this book touches on several very important aspects of computers 
in business—the control features of systems, operations research and management 
science related to system definition, organization and administration of computer 
activity, company planning for computers, etc. There is, however, no originality, 
no profundity, no sureness of touch, no indication of “battle action’? with these 
matters. Instead, there is a hollow reflection of culled magazines. There is, for ex- 
ample, a dearth of case history presented from the author’s vast experience, the 
experience which qualified him to write this book. 

Because there is as yet an extremely limited literature dealing with business 
systems design in general and with principles of the use of computers in such 
systems in particular, this text should prove useful, say, at the level of undergradu- 
ate students in business administration or of the high-speed orientation courses con- 
tained in the “management development programs” of many large companies. 
Better to serve this purpose, there are a sizeable glossary, an extensive collection 
(already badly out of date) of data on commercially available computers after the 
manner of the ONR and BRL surveys, historical material, bibliographic references 
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partitioned by chapter content, and even related homework assignments. The 
popular style should contribute to the book’s palatability as an undergraduate text. 


H. N. Lapen 
The Chesapeake and Ohio Railway Company, 
Terminal Tower, Cleveland 1, Ohio 


46|Z|.—Lovuis CourFricnaL, Les Notions de Base, Gauthier-Villars, Paris, 1958, 60 
p., 24 em. Price $1.55. 


The title of this little pamphlet needs to be interpreted. At the top of the title 
page is Information et Cybernetique, Collection Internationale. 

On a fly leaf are the names of three projected items of the Collection: Boulanger, 
L’ Automation; Prudhomme, Construction des machines automatiques; and Ducasse, 
L’expression des connaissances techniques. 

In the light of these the title of this pamphlet becomes clear, namely, “the basic 
notions of information theory and cybernetics.” 

Monsieur Couffignal merely lists and comments on the basic notions; there is 
nothing original. In fact, a fair part of the text consists of quotations from M. 
Couffignal’s other publications. 

Wiener first published the doctrine of cybernetics in 1948. The term comes from 
the Greek word for pilot. Aboard ship the captain decides where to go and the 
helmsman works the mechanism of the ship, but the pilot between them contrives 
that the work of the helmsman shall achieve the object of the captain. Cybernetics 
is the study of means—principally feed-back—of reaching prescribed goals. Subse- 
quently Wiener published The Human Use of Human Beings about the social im- 
plications of cybernetics. 

Couffignal contests the scientific character of these two books. It is in the 
analysis of human action that Couffignal hopes to find essential principles. Man lives 
in an environment. His actions have an object, to affect the environment. Each 
action is preceded by preparation for the action, a program. Preceding the program 
is a decision to act, based on judgements of values. Thus a definite structure ap- 
pears. 

Having acted, there are three possibilities. First and simplest, the environment 
is affected as foreseen in the program. Or, second, the environment reacts in an un- 
foreseeable fashion, but according to known laws. A simple example is the tempera- 
ture regulation of a building; because of unknowns such as the number of times 
doors will open or the number and activity of the occupants, it is not possible to 
prescribe directly how to regulate the furnace. But a simple feed-back makes pos- 
sible a program that accomplishes the desired objective. Third, the environment 
reacts according to unknown laws in an unforeseeable way, as a pursued butterfly. 

Each of the situations above requires information about the environment; the 
less the reaction of the environment is understood, the more information is needed 
to build an effective program. Thus information is seen to be an essential of cyber- 
netics. 

The word “information” comes from the French, and originally designated the 
action of giving a form; it still has this connotation in its legal usage. The concept 
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“quantity of information” is new. A specific occurrence of information has a phys- 
ical form which is irrelevant to its meaning. For instance, a television program is 
fed over a microwave network as modulation of super-high-frequency radio waves, 
at a transmitting station it may be remodulated onto ultra-high-frequency or very- 
high-frequency waves, or recorded on magnetic tape. At a receiver it is viewed as a 
picture on a kinescope and it may be photographed. Each of these physical forms is 
different and is called a “support” by Couffignal. The common property of all the 
possible forms of an information he calls its ‘“‘semantic.” 

Couffignal argues that cybernetics is not a science. Neither is it a technology 
nor an applied science. He concludes that it is an art, the art of insuring effective- 
ness of action. The material with which it works is information. 

He conceives for each technique a set of systems in which the technique is effec- 
tive. This set of systems is called the “domain of effectiveness” of the technique. 
It is intended that subsequent works in the series of which this is the first will 
report progress in exploring these domains. Some possibilities are: “subjective man”, 
who reacts to his concept of his environment instead of the environment itself; 
“social man’’, who reacts to a social environment, or who reacts en mass as a social 
unit. The theory of governors may yet apply to Governors with a capital G. There 
is the domain of machines, especially information machines. A parallel domain is 
that of automation, the replacing of human beings by machines. There is the 
domain of models and simulation, which includes mathematical models (such as 
the Maxwell equations) of the physical world. And then there is the domain of 
knowledge. This is clearly basic. In fact, it is overwhelming that Couffignal stakes 
out such a broad claim. 

The last chapter is a bald list of concepts. It is in six categories: human beings; 
human actions; information; mechanisms; analogues; afid mentality. 

Evidently M. Couffignal believes that he is starting something that will have 
far-reaching consequences. He has mapped out a broad outline that will take a 
generation to fill in. But he has only outlined broad areas with vague boundaries; 
subsequent items in this collection will need to be much more precise if they are to 
have scientific value. 


H. CAMPAIGNE 
National Security Agency, 
Ft. George G. Meade, Maryland 


47|Z|.—J. von NeuMANN, The Computer and the Brain, Yale University Press, New 
Haven, Conn., 1958, xiv + 82 p., 21 cm. Price $3.00. 


This small volume constitutes the last contribution of one of the great scientists 
of our time, John von Neumann. The material was prepared for delivery at Yale 
University during the spring of 1956 in the Silliman Lectures series. However, the 
lectures were never given, owing to the increasing severity of the illness which 
finally took von Neumann’s life on February 8, 1957. The subject is one which 
attracted his interest for a number of years. It is a subject in which he was eminently 
qualified, by virtue of his great genius and his contributions in the fields of high- 
speed calculators and the logic of automata. In spite of the preliminary nature of 
this work, it is destined to become the nucleus of a new field of research which will 
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challenge the minds of men for many years to come—the comparative study of the 
human brain and man-made automata. 

The book is divided in two parts. In the first part von Neumann discusses the 
basic principles underlying the design of modern computing machines, both analog 
and digital. Some of the properties of analog computers receiving special attention 
are: (1) their continuous but approximate method for representing quantitative 
information; (2) their capacity to obtain valid solutions to many classes of mathe- 
matical problems in spite of their limited accuracy, which is, at best, 0.01%; and 
(3) the use of basic operations more advanced than the four arithmetic operations 
in some analog calculators such as the “Stieltjes” integral in the differential analyzer. 
The characteristics of modern electronic digital calculators are discussed in some- 
what greater detail. Of special interest for comparison with the human brain are: 
(1) the quantized or “marker”? method for representation of information; (2) the 
conventional basic arithmetic and logical operations used; (3) the organization of 
the controls of digital calculators, which is governed by strict logical rules and 
prescribed sequences; (4) the high precision available in these computers, which 
is approximately 10-" in large-scale computers currently in operation; (5) the 
organization and capacity of the memory components; (6) the great arithmetic and 
logical “depth” required for the solution of most problems on digital computers 
(i.e. the very large number of operations used in sequence); (7) the use of the 
memory for storing instructions as well as other information, resulting in the cap- 
ability of digital computers to modify or operate on instructions; and (8) the use of 
interpretive subroutines or ‘“‘short codes” for carrying out functions more advanced 
than the basic operations built into the hardware of a computer. One of the key 
points made by von Neumann is that the very high precision requirement of digital 
computers is dictated by the inherent deterioration of accuracy resulting from the 
very large number, or great “depth,” of the mathematical and logical operations 
used in sequence. 

In the second part von Neumann compares the functioning of the human brain 
with the operation of a modern computer, bringing out the areas of “similarity and 
dissimilarity between these two kinds of automata.” He begins by describing the 
known properties of a neuron and the apparent digital character of its mechanism 
for receiving and transmitting pulses. Some of the more superficial comparisons be- 
tween the human brain and its man-made counterpart are: (1) speed—there is a 
factor of about 10* to 105 in favor of the man-made components; (2) size—there is a 
factor of about 108 to 10° in favor of the human brain (the number of neurons esti- 
mated in the central nervous system is of the order of 10'°); (3) energy disstpation— 
there is a factor of 10° to 10° in favor of the human components; and (4) accuracy— 
there is a factor about 108 in favor of the artificial componentry. 

Looking deeper into the more intrinsic areas of comparison, von Neumann is 
led to the conclusion that the basic internal language used by the brain is un- 
doubtedly quite different from the mathematical language with which we are 
acquainted. He arrives at this conclusion primarily on the basis of the argument 
that the information stored in the brain lacks sufficient accuracy to enable it to 
carry out mathematical and logical processes in such “depth” as would be required 
if the language used were based on conventional mathematical symbols. Von Neu- 
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mann conjectures that the language used by the brain is probably statistical in 
nature in which correlation processes play an important role. (See The Perceptron— 
A theory of statistical separability in cognitive systems, Cornell Aero. Lab. Inc., 
Report Nos. VG 1196 G-1, 1958 and 1196 G-2, by F. Rosenblatt.) He concludes the 
book with the remark: “‘Thus logics and mathematics in the central nervous system, 
when received as languages, must structurally be essentially different from those 
languages to which our common experience refers. 

“Tt also ought to be noted that the language here involved may well correspond 
to a short code in the sense described earlier, rather than to a complete code: when 
we talk mathematics, we may be discussing a secondary language, built on the 
primary language truly used by the central nervous system. Thus the outward forms 
of our mathematics are noi absolutely relevant from the point of view of evaluating 
what the mathematical or logical language truly used by the central nervous system 
is. However, the above remarks about reliability and logical and arithmetical depth 
prove that whatever the system is, it cannot fail to differ considerably from what we 
consciously and explicitly consider as mathematics.” 

(Courtesy of Applied Mechanics Reviews) H.P. 


48[Z)._Cuares V. L. Smiru, Electronic Digital Computers, McGraw-Hill, New 
York, 1959, xi + 443 p., 24 em. Price $12.00. 


In the preface the author states: “This is not a treatise on digital computer 
engineering, nor on the other hand an exhaustive treatise on logical design. The 
reader will find considerable discussion of circuits and components—enough, I 
hope, to give him a reasonably complete understanding of various ways in which 
the usual functions of a computer, such as memory, control, the performance of 
arithmetic, and input and output, can be realized physically. But I have not at- 
tempted a treatment sufficiently detailed to provide design information. The reader 
will also find sufficient information on computer arithmetic and instruction codes to 
provide him with a basic understanding of these matters, but here again I have not 
attempted the detailed treatment that a logical designer would demand, and I have 
for brevity considered only machines using binary arithmetic. My purpose has been 
to provide the reader with sufficient information to understand how digital com- 
puters function.” 

No book of this size could possibly cover every phase of computers; however, 
for the numerous topics chosen, the author’s objective has been met. Mathema- 
ticlans, programmers, and engineers who have limited knowledge of the basic 
principles should find this volume especially useful in extending their understand- 
ings. As stated previously, this is not a treatise, and its organization is such that 
it is probably most useful as a reference text. There are many additional references 
footnoted throughout the book for those desiring additional details. The usefulness 
of the subject index, however, suffers somewhat because of its brevity. The bulk of 
the material has been written about computers and techniques of 1956 and earlier, 
but this need not detract in any way from the value of the material. 

Chapter titles are: 

1. Digital-computer Arithmetic (Number Systems and Their Machine Repre- 
sentation) 

2. Instruction Codes 
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3. Some General Considerations on Systems (Basic Functions and Inner Struc- 
ture) 

4. Basic Logie Circuits and Their Representation (Logic or Switching Circuits) 

5. Static Memory Cells (Circuits and Devices Having Two Stable States) 

6. Dynamic Memory Cells (Circuits and Devices Which Require Recirculation 
or Regeneration) 

7. Higher-order Logic Circuits (Boolean Algebra, Switching Matrices) 

8. Shifting Registers (Basic Storage Combined into Complex Structures) 

9. Counters 

10. Adders and Accumulators 

11. Large-scale Memory Devices I (Magnetic Drums, Ultrasonic and Mag- 
netostrictive Delay Lines) 

12. Large-scale Memory Devices II (Electrostatic Tubes and Magnetic Cores) 

13. The Arithmetic Unit 

14. The Memory 

15. The Central Control 

16. Input and Output 

17. Superspeed Computers (NBS, LARC, STRETCH, ILLIAC II). 
Gorpon D. GoLpsTEIN 


Office of Naval Research, 
Washington, District of Columbia 








TABLE ERRATA 


272.—H. M. Navtican Atmanac Orrice, Planetary Coordinates, 1960-1980, 

H. M. Stationery Office, London, 1958. 

P. 145, line 18; for n’ read n?. 

P. 146, line 14; for 59.8 read 58.1. 

This erroneous value for m is also on page xiv; the correct value brings the 
results of the example into still closer agreement with those on page xiii. 

J. G. PorTER 

H. M. Nautical Almanac Office, 


Royal Greenwich Observatory, 
Hailsham, Sussex, England 
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NOTES 


Vychislitel’naia Matematika 


The new Russian Journal Vychislitel’naia Matematika has recently appeared. 
The subjects treated in the first two issues are in the fields of Numerical Methods 
and Computing Machines (both digital and analogue). We reprint below the table 
of contents of the first two volumes. 

1. Izdatel’stvo Akademii Nauk SSSR, Vychislitel’nyi Tsentr, Sbornik 1, Moskva, 
1957. 

L. A. Liusternik; The finite-difference analog of Green’s function in the three- 
dimensional case. 

iv. V. Vorob’ev; The method of moments in the problem of the vibrations of 
linear systems. 

E. A. Volkov; Investigation of a method of improving the accuracy of the method 
of nets for solving Poisson’s equation. 

V. K. Saul’ev; On a class of elliptic equations, solved by use of the method of 
finite differences. 

V. K. Saul’ev; On an estimate of the error in getting characteristic functions by 
the method of finite differences. 

A. I. Vzorova; On the construction of polynomials orthogonal on a family of 
ellipses. 

[A. I. Alikhashkin; A method of calculating the discharge for forced flows through 
an imperfect slit. 

fA. I. Alikhashkin; Solution of the problem of the imperfect slit by the method 
of straight lines. 

G. S. Khovanskii; On the replacement of the logarithm function by a power in 
making approximate nomograms. 

B. P. Kapitsa; Mechanical computation of the harmonic adjoint function. 

2. Izdatel’stvo Akademii Nauk SSSR, Vychislitel’nyi Tsentr, Sbornik 2, Moskva, 
1957. 

M. R. Shura-Bura; Approximating functions of many variables by functions 
each of which depends on one variable. 

P. I. Chushkin; The flow around ellipses and ellipsoids in sonic gas flow. 

O. N. Katskova and [U. D. Shmyglevskii; Axisymmetric supersonic flow of a 
freely expanding gas with a plane transitional surface. 

V.S. Linskii; Computation of elementary function on automatic computing ma- 
chines. 

M. G. Rappoport; Computation of finite differences on punched-card machines. 

B. M. Drozdov and M. G. Rappoport; Coding of operations on the electronic 
computer EV80-3. 

A. G. Gesse; Electrical resonance networks for solving systems of linear equations. 

G. 8S. Khovanskii; Methodology of constructing nomograms with a triangular 
(hexagonal) transparency. 

E.I. 





































NOTES 


New Journals 


The publication of two new German journals in the fields of numerical analysis 
and data-processing has been announced recently. The journal in numerical analysis 
has the title Numerische Mathematik and is published by Lange & Springer, 
West-Berlin. The editors are: R. Sauer, E. Stiefel, J. Todd, and A. Walther. The 
first issue contains the following papers: 

1. Uber diskrete und lineare Tschebyscheff-Approximationen. 

2. On certain methods for expanding the characteristic polynomial. 

3. Orthogonal polynomials in several variables. 

4. Report on the algorithmic language ALGOL. 

The new journal on data-processing is published by Friedr. Vieweg & Sohn,’ 
Braunschweig under the editorship of H. K. Schuff. Its title is Elektronische Daten- 
verarbeitung. The first issue contains a number of papers on programming and ap- 
plication of digital computers with particular emphasis on their use in the manage- 
ment data analysis field. 


H.P. 








